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SUMMARY 


A computational study has been undertaken to study the performance of advanced 
phenomenological turbulence models coded in a modular form to describe incompressible turbulent 
flow behavior in two dimensional/axisymmetric and three dimensional complex geometry. The 
models include a variety of two equation models (single and multi-scale k -e models with different 
near wall treatments) and second moment algebraic and full Reynolds stress closure models. These 
models were systematically assessed to evaluate their performance in complex flows with rotation, 
curvature and separation. The models are coded as self contained modules that can be interfaced 
with a number of flow solvers. These modules are stand alone satellite programs that come with 
their own formulation, finite-volume discretization scheme, solver and boundary condition 
implementation. They will take as input (from any generic Navier-Stokes solver) the velocity field, 
grid (structured H-type grid) and computational domain specification (boundary conditions), and 
will deliver, depending on the model used, turbulent viscosity, or the components of the Reynolds 
stress tensor u t uj . There are separate 2D/axisymmetric and/or 3D decks for each module 
considered. 

The modules are tested using Rocketdyn's proprietary code REACT. The code utilizes an efficient 
solution procedure to solve Navier-Stokes equations in a non-orthogonal body-fitted coordinate 
system. The differential equations are discretized over a finite-volume grid using a non-staggered 
variable arrangement and an efficient solution procedure based on the SIMPLE algorithm for the 
velocity-pressure coupling is used. The modules developed have been interfaced and tested using 
finite-volume, pressure-correction CFD solvers which are widely used in the CFD community. 
Other solvers can also be used to test these modules since they are independently structured with 
their own discretization scheme and solver methodology. Many of these modules have been 
independently tested by Professor C.P. Chen and his group at the University of Alabama at 
Huntsville (UAH) by interfacing them with own flow solver (MAST). 
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CHAPTER 1 


Introduction 


1.1 Background 

Computational Fluid Dynamics (CFD) has been used extensively for the last decade or so in 
analyzing complex flow phenomenon for many industrial applications, such as combustion and 
turbomachinery. Most flows of practical interest are turbulent and for many of them, relatively 
simple prediction methods are sufficient to produce results of engineering accuracy. For others, 
mainly flows in complex geometry with large body forces such as curvature, rotation and 
separation, more complex prediction methods are required. 

With advancing state-of-the-art of computer technology, the range, size and complexity of flow 
models being applied have increased. Users have become more sophisticated and there is a 
constant demand for improvement. CFD codes have adapted to this demand and many general- 
purpose computer codes have been developed and used. As these general purpose codes become 
larger, their code structure becomes sophisticated and in general this structure can be divided into 
three main areas; 

1) Numerical algorithms which include discretization methods and solution techniques. 

2) Methods of dealing with complex geometry, such as grid generation, structured or 
unstructured grids. 

3) Physical models which include turbulence models, porosity, combustion kinetics, 
multi-phase flows, etc. 

It seems, therefore, that the practicing engineer must have the knowledge of all these elements of 
the CFD program in order to successfully utilize the code. Modularization of the code structure 
may then become necessary in order to obtain the maximum benefits from these general-purpose 
CFD codes. This means developing individual modular routines for the solver and other physical 
models. If such modules are successful they would allow users to concentrate their talents on 
developing and improving physical hypothesis such as turbulence models that can be easily tested 
using these modules. 
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In general, the physics of turbulence can be captured by solving the full time-dependent Navier- 
Stokes equations in what is termed as Direct Numerical Simulation (DNS). However, DNS is not 
practical for engineering purposes mainly because it is restricted to flows at low Reynolds 
numbers. Large Eddy Simulations (LES) are now competitive with DNS in accuracy at an order of 
magnitude less cost, however, it is still expensive for routine engineering calculations. Therefore, 
current engineering prediction methods are based on Reynolds-averaged equations, with models 
for the unknown Reynolds stresses which appear as the result of time-averaging the nonlinear 
Navier-Stokes equations. These models fall mainly into three categories; "eddy-viscosity" models, 
where a relation between the Reynolds stresses and mean velocity gradients at the same point in 
space is sought. Algebraic stress models, where the Reynolds stresses are expressed as an 
algebraic relation of turbulence production and dissipation. Reynolds stress models where the exact 
partial differential equations for the Reynolds stresses are solved after closing the higher order 
terms. These transport equations account for the dependence of Reynolds stresses on the history of 
the flow and should perform better than the eddy-viscosity models. 


1.2 Outline of the Present Study 

In the present work, phenomenological, single-point turbulence models coded in a modular format 
are developed as self-contained code decks that can be interfaced with a number of flow solvers to 
analyze turbulent flows in complex 2D/axisymmetric or 3D geometry. These modules are validated 
using Rocketdyn's REACT code and are independently tested at UAH using own code MAST. 

The models that are developed in a modular form include; 

1 _ 2D/axisymmetric single-scale k-£ model with three options for near wall treatment that include, 

- Standard Launder and Splading wall functions. 

- Chen and Patel two-layer model. 

- Lam and Bremhorst low-Reynolds number model. 

2. 2D/axisymmetric multi-scale k-e model with the standard wall function and Chen & Patel two- 
layer near wall treatment. 

3 . 2D/axisymmetric implicit algebraic stress model (ASM) based on the original work of Rodi. 

4. 2D/axisymmetric full Reynolds stress turbulence model (RSM) based on the simplified linear 
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second moment closure model of Launder, Reece and Rodi (LRR) second moment closure. 


5 . 3D standard k-e turbulence model with wall function and two-layer near wall treatments. 

6. 3D algebraic stress model (ASM). 

Each model is coded as a self contained, stand alone module deck that can be interfaced with a 
number of CFD solvers to analyze turbulent flows in complex geometry. The user can use these 
modules without concern as to how they are implemented and solved. The input to the modules are 
the mean flow variables, boundary and geometric information which are to be provided by a mean 
flow solver. The output of the module are the turbulent eddy-viscosity for the eddy-viscosity 
models and the Reynolds stresses for the second moment closure models. Moreover, source terms 
which are needed for the mean flow calculations are calculated and must be passed to the main 
solver. The modules are tested using the finite-volume REACT code and the results compared with 
available experimental data. 

Full details of each module are given in the next chapters. Chapter 2 discusses the theory and 
model equations for the two-equation k-e model used in the 2D/axisymmetric module deck. The 
module is evaluated with a number of benchmark problems and detailed description of the module 
variable names together with the input/output structure are given in appendix A. The complete 
listing of the module is provided at the end of the chapter. Similarly, chapter 3 discusses the theory 
and model equations for the 2D/axisymmetric multi-time-scale k-e model. The 2D/axisymmetric 
Algebraic stress module is presented in chapter 4 and chapter 5 discusses the 2D/axisymmetric 
Reynolds stress module deck. Full description of the 3D k-e turbulence model is given in chapter 
6 and chapter 7 presents a full description of the 3D algebraic stress model together with module 
description and code listing in the appendix. Finally in chapter 8, copies of related turbulence work 
that are presented or published elsewhere are attached. 
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2.1 Introduction 


In this section a description of the standard k-e turbulence model that is coded as a self contained 
computer program to compute turbulent flow quantities in two-dimensional planar or axisymmetric 
geometry is given. Detailed description of the module structure, variables used and how to 
interface the module with CFD flow solvers are given in Appendix A. The module has been tested 
as a separate self-contained unit using the REACT code [1] and was independently tested at the 
University of Alabama at Huntsville (UAH) using own code (MAST). 

2.2 Theory and Model Equations 

The k-e turbulence module is based on the widely used single-scale two equation k-e turbulence 
model (Jt is the turbulent kinetic energy and e is the energy dissipation rate). The model developed 
originally by Launder and Spalding [2] was successful in providing good predictions for a wide 
range of turbulent flows. The k and e -equations can be derived from the transport equations for 

the Reynolds stresses assuming fully turbulent flow. 

For low-Reynolds number flows close to solid boundaries, adjustments to the model are needed to 
bridge the viscous dominated sublayer region with the fully turbulent flow region. The success of 
the wall function method depends on the universality of the turbulent flow structure near the wall. 
In many complex flows, however, the flow field near the wall has to be determined accurately and 
the traditional wall-function method is not satisfactory. This is because the specification of all 
turbulence quantities in terms of the friction velocity fail at separation where the flow near the wall 
is no longer controlled by the wall shear stress. Patel et al [3] assessed the relative performance of 
various models which describe the near-wall flows and found that there are still areas of 
improvements needed to accurately model flow behavior near the wall. 

Jones and Launder [4] extended the original k-e model to the low-Reynolds number form which 
allowed the calculation to be performed all the way to the wall. Numerical difficulties of accurately 
resolving the large gradients close to the wall necessitates resolving the wall region with a very fine 
grid structure. Chen and Patel [5] introduced a method to resolve the near-wall region which 
combines the standard k-e model with the one-equation model of Wolfshtein [6] near the wall. In 
this "two-layer" model an algebraically prescribed eddy-viscosity for the wall region is coupled to 
the k-e model to describe the details of the flow in the vicinity of the wall. 
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Momentum and continuity equations are solved up to the wall and this reduces the physical 
uncertainties of near-wall turbulence and the numerical difficulties of resolving the very large 
gradients of turbulence parameters. 


For an incompressible, steady and axisymmetric turbulent flow, the Reynolds averaged momentum 
and continuity equations can be expressed in a generalized form as; 


d(puO) 

Bx 


+ 7 J: (pvr&) = J 7 (r& x ^) +7 ( rF ®r^) + s * 


( 1 ) 


where <Z> is the dependent variable, which stands for = u t v, w for the axial, radial and 
tangential velocities respectively, p is the fluid density, JHcj^ and r<£> r are exchange coefficients in 
x and r -directions, respectively, and S<p is the source term for the variable 0. 


The source terms for the dependent variable are: 

• Axial direction, 0-u, r<p x = 2/i ef IcEy = ji e and 

o BP 1 B (llr Bv ) 

Su = + TdF ( ^ er Bx } 
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where p e is the eddy viscosity and P is the pressure 


• Radial direction, <Z>=v, r<p x = p e , r& r = 2p e and 
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Equations 2, 3, and 4 above are the momentum equations that are solved by the CFD solvers. 
However, in order to close the equations and determine the eddy viscosity different turbulence 
models are used. 



The present module utilizes the k-e model. In this model two equations for the turbulent kinetic 
energy k and its dissipation £ which have the same general form as equation ( 1 ) are solved. 

For the turbulent kinetic energy equation 

p t 

0 = k, r<p = r<p r = p + and S& = G - pe ( 5 ) 

°k 


For the energy dissipation equation 


O = £ r<p x = r&= p+ — and S& = £ (Cjf jG - Cfa P £ ) 

°£ 


( 6 ) 


where 0 £ and CTg are turbulent Prandtl/Schmidt numbers for k and e respectively, and G denotes 
the rate of production of the turbulent kinetic energy and is expressed as: 


G - p e { 2 [(^) 2 + (^r) 2 + 


dx } 


[ dr > 




(7) 


where fi is the dynamic viscosity, and p t is the turbulent viscosity, 

Pt = Ctfn p — (8) 

£ 

and p e = p + p t 

C^, C], C 2 , Ok and a e are constants whose values are 0.09, 1.44, 1.92, 1.0, 1.0, 
respectively and/ 7,/2 and/ m are damping functions. 

Near a wall, turbulent flow can be divided into two regions, the inner viscous sublayer where low 
turbulent Reynolds number effects are important and the velocities decrease rapidly to zero at the 
wall, and the outer fully turbulent region. The successful application of the k-e turbulence model 
for many complex flows depends to a large extent on how accurately the flow field near the wall is 
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determined. In the present module three different models are used to treat this thin sublayer region, 
they include; 

Wall function method, where 

u + =y + at y + < 11.6 (9) 

u + = — In ( Ey + ) at y + > 11.6 (10) 

K 

where, u + = — , y + = and u T = V z 

u x y 

z w is the wall shear stress which can be determined from 


Tw = for y + ^ 11.6 
o 


( 11 ) 




K p Ujy fc 0 ,5 

In [E C /- 25 p 8k 0 5 /fi] 


for y + > 11.6 


( 12 ) 


Here, Up denotes the velocity component parallel to the wall at the first grid point p from the wall. 

8 is the normal distance from the wall and K is a constant = 0.42. 

In this approach, k and £ equations are solved with fp =fj -f2 - 1, only in the fully turbulent 
region beyond some distance from the wall. Boundary conditions i.e., velocity components and 
turbulent parameters at that distance are specified in terms of the friction velocity u x . 

In the low-Reynolds number model, the flow is resolved all the way to the wall with a very fine 
mesh. Many models have been proposed that are based on the k-e model and differ mainly in the 
choice of the damping functions fp , fj and /2 to bridge the gap between the sublayer and the fully 
turbulent region. The model due to Lam & Bremhorst [7] is used in this work, where; 

fp = [ 1 -exp (-0.016 R y )] 1/2 (1+^-) 
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fl = 1 + ( - j jp) 3 and f 2 = 1 - exp ( - R 2 { ) 


k ,/2 y 

where, R v = ^ and Rf 

J v 


k 2 

— are turbulent Reynolds number, 

v e 


These damping functions tend to unity with increasing distance from the wall. 

In the two-layer model due to Chen and Patel [5], a simple algebraically prescribed eddy-viscosity 
model for the wall region is coupled to the k-e model for the outer flow to describe the flow 
details. Unlike the low-Reynolds number model that requires the solution of transport equations 
for both it and e all the way to the wall, the one-equation model requires the solution of only the 
turbulent kinetic energy equation in the sublayer region while algebraically specifying the eddy 
viscosity and energy dissipation. 

i A/2 iJ/2 

v t = CpL-— and £ = — . 

L H l £ 

The length scales and L £ contain the necessary damping effects in the near-wall region in 
terms of the turbulence Reynolds number Ry. 

L y. = c iy [1 -exp(-R y /Afi)] ( 13 ) 

L £ = Cjy [1 - exp (-Ry/A^] (14) 

and L £ become linear and approach C } y with increasing distance from the wall. 

Cj = kC°' 75 and A s = 2Cj. Chen and Patel [5] used A^ = 70. 

The damping effects decay rapidly with distance from the wall independent of the magnitude of the 
wall shear stress. The matching between the one-equation and the standard k-£ models is earned 
along prescribed grid lines where Ry ~200. 

For flows in rotating ducts a modification was made by Chen and Guo [8] to reflect the effects of a 
system rotation on the length scales andL £ , as; 
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Ljj = Ljjq [ 1.0 + 1.3(0.4^-0.8Q)Q(^) 2 ] 15 
L e = L eo [1.0+1.3(0.4^-0.8Q)£2(^) 2 ] 0 - 5 
Moreover, the function /2 in the dissipation equation is modified to 

/2 =f2 + M 

where Ri is a Richardson number to reflect the effects of streamline curvature due to rotation and 
is defined as 

R t = ( 0.4 (O k - 0.8 Q k ) Q k (— ) 2 

£ 

dU- 

where 0) k = e ij k is the local mean vorticity . 

The above modification to account for streamline curvature and rotation seemed adequate in the 
framework of two equation k-e modeling. Other modifications have also been considered but not 
implemented in this module and can be referred to in Hadid and Sindir [9]. 


2.3 Module Evaluation 

The single scale k-e turbulent module was evaluated by comparison with published experimental 
data. One of the test problems considered is the two dimensional incompressible turbulent flow 
over a backward facing step with and without rotation (see figure 1) to compare with the 
experiment of Rothe and Johnston[10]. While the mean flow is in the x -y plane, the channel is 
rotated with constant angular velocity Q about the z -axis. The ratio of the channel width to the 
step height is very large so that the secondary flow can be ignored, which made the flow remain 
two dimensional. The channel height to step ratio was set to 2 and the inlet channel height ( h ) 
equals to the step height ( H ). The Reynolds number based on the uniform inlet velocity was about 
5500. The rotation number ( Ro = Qh/U ) was varied between +0.06 and -0.06. 

The streamline patterns for the three different rotation numbers Ro = -0.06, 0.0, +0.06 by using 
the three different wall treatments are shown in figures 2-4. In each figure, the upper (a) and lower 
(c) parts correspond to Ro = +0.06 and Ro= - 0.06 respectively. While the middle part (b) is the 
non-rotating case. It is observed that the streamline patterns are influenced by the system rotation. 
Suction side step extends the recirculation zone and the pressure side step reduces the recirculation 
zone. The reattachment length for Ro= - 0.06 using the wall functions is larger compared to the 
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experimental results. This is due to the fact that no Coriolis effect is accounted for in the law of the 
wall. The predicted variation of reattachment length with Ro (figure 5) shows reasonable 
correlation with the experimental data of Rothe and Johnston [10]. 

The single scale k-e model using three different wall treatments with rotational stress generation 
terms embodied seems to capture the main effects of system rotation on turbulence structure, i.e. 
the suppression of turbulence level with clockwise rotation and enhancement of turbulence level 
with counterclockwise rotation The effects are also noticeable in the corresponding increase in the 
reattachment length with clockwise rotation and its decrease with counterclockwise rotation. 

The other two test cases were those of Daily and Nece [11] where rotating disk cavity circulation 
and secondary flows are induced by a rotating wall, and Roback and Johnson [12] for a confined 
double concentric jets with a sudden expansion. Flow swirl in this case is induced by imposing a 
tangential velocity component at the outer jet. Figure (6) shows the two-dimensional axisymmetric 
rotating lid cavity of Daily and Nece. The flow is bounded by a disk (rotor) and a stationary end 
wall (stator) of a chamber. The ratio of the axial clearance between the rotor and the stator (s ) to 
the radius of the disk (a ) is 0.0255. The disk rotates with a rotational Reynolds number 
R=4.4xl0 6 defined as R = Oa 2 /v, where Q is the disk rotational speed and v is the kinematic 

viscosity. 

Computations were performed on a 33x75 grid with different grid clustering near the walls for the 
different near-wall models. Figure (7) shows the velocity vectors at the top region of the cavity 
using the wall function model. Centrifugal forces move the fluid radially outward on the disk, 
axially away from the disk on the wall casing, and radially inwards on the stationary end wall. 
Figure 8, shows the axial variations of the radial velocity component at a radial position r/a=0.765. 
The agreement is fair with some discrepancy for all near-wall models close to the rotating disk. 
Figure (9), shows the axial variation of the tangential velocity component at the radial position 
r/a-0.765. At the rotating disk (x=0), the tangential velocity approaches the value (aO ). The two- 
layer near wall model seem to offer closer agreement with the data than the other two models. The 
presence of corner regions presents a difficulty in defining the normal distances used in the 
definition of turbulent Reynolds number. In the present analysis, values of the normal distance 
were based on the normal distance to the nearest solid boundary. 

Predictions of the experiments of Roback and Johnson [12] have been presented by several 
workers, e.g. Sloan et al. [13] and Durst and Wennerberg [14], Unfortunately, inlet flow profiles 
were not provided in the experiment. Therefore, the present calculations were started at the 
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expansion plane using the measured velocity profile at 5 mm downstream of the expansion after 
some adjustments near the edges of the coaxial jets. Measurements of main turbulent intensities 
were used to calculate inlet values of the turbulent kinetic energy. Energy dissipation rate was 
estimated from e = C^k 3/2 /L, where L is a length scale of turbulence at the inlet of the order of 

10' 4 m. 

Figure 10, shows an illustration of the test chamber geometry. The chamber diameter is about 
twice the secondary tube diameter. The exit from the 8-bladed, 30°, free vortex swirl generator is 
located approximately 0.005 m upstream from the confluence plane. 

A prominent phenomenon in axisymmetric swirling flows in such geometry is the "bubble” or 
vortex breakdown which has been studied extensively [15-18]. In the present numerical simulation 
of the experiment, a 150x100 grid nodes was used with different clustering on the walls for the 
different near-wall models used. Figure 11, shows the velocity vectors indicating the presence of a 
closed recirculation zone at the center with additional zones at the comer downstream and between 
the inner jet and the outward diverted secondary jet. The figure also shows flow diversion 
outwards with high gradients characterized by large turbulent shear and fluctuation levels. 
Comparisons were made of the radial variations of flow variables at two axial locations, x =0.02 5m 
upstream of the vortex bubble and x=0.102m located inside the bubble. Figure (12), shows the 
radial variation of the axial velocity profile at x=0.025m using the wall function, two-layer and the 
low Reynolds number models. Fair agreement by the different models is shown. They also seem 
to predict small negative velocities at a radial position r~0.0153m (the interface between the inner 
and outer jets), slightly under predicted in strength and width. Figure (13) shows the radial 
variation of the axial velocity profiles at x=0.102m. The two-layer model shows a better agreement 
with the experimental data. 

Radial variations of the tangential velocity at x=0.025m is shown in figures 14. The figure shows 
that the two-layer model offers better agreement with the experiment as compared with the wall 
function or the low Reynolds number models. 

In general, the calculations shown above indicate that the two layer model seem to offer a better 
comparisons with the experimental results. The three near- wall models are built in the standard 
two-dimensional/axisymmetric k-£ turbulence module. The structure of the module will be 
discussed next together with the details of interfacing with a flow solver and descriptions of 
variables. 
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Figure 1. Rotating backward facing step 
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Figure 2. Stream-function contours using 
wall function near wall treatment 
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Figure 3, 


Stream-function contours using the 
two-layer near-wall treatment 
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Figure 4. Stream-function contours using the low- 
Reynolds number near-wall treatment 
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Figure 5. Reattachment length as a function of the 
rotation number 
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Figure 10. Roback & Johnson's swirling coaxial 
jets discharging into an expanded duct 
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Figure 11. Velocity vectors 








APPENDIX A 


2D/Axisymmetric k-e Turbulence Module Deck 
A.l Introduction 

In an attempt to modularize the k-e turbulent physical model -a difficult task as many CFD users 
may know. A self-contained, stand-alone turbulence module has been constructed that computes 
turbulent flow quantities using the standard k-e turbulence model. The module is structured to be 
flexible with options for three near- wall treatments. It can be easily accessed by the user and 
interfaced with own CFD solvers to calculate turbulent flows. 

It is hoped that the program is sufficiently "full proof and user friendly. However, care must be 
exercised to identify the limitations of the module to be compatible with the flow solver. Module 
capabilities and input/output structure is described next in details followed by a FORTRAN listing 
of the module. 


A.2 Program KEMOD 

This is basically the solver for the k and e - transport equations. It reads through its argument list 
different variables from the calling flow solver. These variables are described below where, each 
variable name ends with either an (I) for Integer variable, (R) for Real variable or (L) for Logical 
variable. 

The flow chart of the program is shown in Figure A.l. It shows the main operations performed by 
the code. 


List of Argument Variable Names 

NIMI Number of cell nodes in the I- or ^-coordinate lines, (input from flow solver) 

NJMI Number of cell nodes in the J- or T|-coordinate lines, (input from flow solver) 
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XR 

YR 

UR 

VR 

WR 

TER 

EDR 

URFKR 

URFER 

PRTKR 

PRTER 

GR 


Grid node locations in the x or ^-direction, dimensioned to XR (NX,NY) (input 
from flow solver) 

Grid node locations in the y or rj-direction, dimensioned to YR (NX,NY) (input 
from flow solver) 

Axial or x-direction velocity (u), dimensioned as UR (NX, NY) (input from 
flow solver 

Radial or y-direction velocity (v), also dimensional as VR (NX,NY) (input 
from flow solver) 

Azimuthal velocity (w), dimensional WR (NX,NY) (input from flow solver) 

Turbulence kinetic energy k, dimensioned TER (NX,NY) (calculated in 
KEMOD and returned to flow solver) 

Turbulent energy dissipation rate e, dimensioned EDR (NX,NY) (calculated in 
KEMOD and returned to flow solver) 

Under-relaxation factor for k -equation (input from flow solver) 

Under-relaxation factor for e -equation (input from flow solver) 

Prandtl/Schmidt number for turbulent energy-equation, assumed known (input 
from flow solver) 

Prandtl/Schmidt number for turbulent energy dissipation equation, assumed 
known (input from flow solver) 

= 1.0 if second order upwinding is desired 
= 0.0 if first order upwinding is used 

(input from flow solver. Usually calculation of k and e are not very sensitive 
to the order of upwinding used) 
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FIR 

F2R 

ITERI 

VISCOSR 

VISR 

AKSIL 

LREL 

LAY2L 

C1R 

C2R 

CMUR 

I2LWI 

I2LEI 

J2LSI 


Mass flux variable at cell faces in x- or ^-direction, dimensioned FIR (NX,NY) 
(input from flow solver) 

Mass flux variable at cell faces in y or r| -direction, dimensioned F2R (NX, NY) 
(input from flow solver) 

Iteration number (input from flow solver) 

Dynamic viscosity (input from flow solver) 

Eddy viscosity, dimensioned VISR (NX,NY) (calculated in KEMOD and 
returned to main solver) 

Logical variable for axisymmetric geometry (AKSIL=-TRUE-) or plain 
geometry (AKSIL=-FALSE-) (input from flow solver) 

Logical variable for Lam & Bremhorst's low-Reynolds number model 
(LREL=-TRUE-) or others (LREL=- FALSE-) (input from flow solver) 
Logical variable for Patel's two-layer model if (LAY2L=*TRUE*) or others 
(LAY2L = -FALSE-) (input from flow solver) 

Turbulence model constant, Ci (input from flow solver) 

Turbulence model constant, C 2 (input from flow solver) 

Turbulence model constant, Cjj. (input from flow solver) 

Grid line location from the west wall in the x-direction for the two-layer model 
(input from flow solver) 

Grid line location from the east wall in the x-direction for the two-layer model 
(input from flow solver) 

Grid line location from the south wall in the y-direction for the two-layer model 
(input from flow solver) 
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J2LNI Grid line location from the north wall in the y-direction for the two-layer model 

(input from flow solver) 

JTBEI Boundary condition flag along east boundary must have one for each boundary 

node set to: 1 -inlet, 2-outlet, 3-symmetry and 4-wall e.g., for an outlet 

boundary condition on the east boundary set JTBEI to NJ*2, and similarly for 
other boundaries, dimensioned JTBEI (NY) (input from flow solver) 

JTBWI Boundary condition flag along west boundary, dimensioned JTBWI (NY) 

(input from flow solver) 

ITBNI Boundary condition flag along north boundary, dimensioned ITBNI (NX) 

(input from flow solver) 

ITBSI Boundary condition flag along south boundary, dimensioned ITBSI (NY) 

(input from flow solver) 

Program KEMOD is interfaced with the main flow solver by a call to KEMOD with its arguments. 
For iterative flow solvers KEMOD is called within the iteration sequence after the solution of the 
momentum equations where the mean velocities are passed to KEMOD. There are different flow 
solvers utilizing different schemes from staggered to nonstaggered grid arrangement and for 
nonorthogonal coordinate system there are at least three alternatives to the choice of the velocity 
components 

i. Cartesian velocity components 

ii. Contravariant velocity components 

iii. Covariant velocity components 

The Cartesian velocity components are the most widely used and have the advantage of simple 
formulation of the governing equations. Whatever the arrangement used, mass fluxes at cell faces 
are required and passed to KEMOD as FIR and F2R in both directions. The location of other 
variables such as k and £ are at the cell center or cell nodes. 
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The module starts by reassigning variable names passed to it from flow solver to names that are 
shared with the different subroutines of the module in a common statement file 
"KEMOD-COMMON Then a check is made if it is the first iteration in which case the grid file 
"GRIDF" is called -after passing the grid node locations XR & YR in KEMOD- in order to 
calculate grid related quantities which will be explained later. The need to call GRIDG can be 
waived if all the grid data are passed to the module. That is all the information about the grid such 
as interpolation factors FX and FY, cell areas (ARE) and volumes (VOL) and normal distances of 
first grid point from grid boundaries (DNS from south boundary, DNN - from north boundary, 
DNW - from west boundary and DNE - from east boundary). 

After this a call to subroutine CALCE is made to calculate the turbulent kinetic energy k (with the 
identifier IPHI=1) followed by a check if the low-Reynolds number model or the two-layer model 
are to be used in which case subroutine TWOLAY is called. The energy dissipation equation is 
solved next by a call to subroutine CALCE again with the identifier IPHI=2. The turbulent 
viscosity is updated next by calling subroutine MOD VIS. A brief description of each subroutine is 
given next. 


A.3 Subroutines 
GRIDG 


Before calling this subroutine, the coordinates of all grid nodes, defined in reference to a fixed 
Cartesian coordinate frame are read. Figure A.2 shows the position of cell and grid nodes. 

This subroutine is called only once to calculate coordinates of grid nodes (intersection of grid lines) 
and geometrical properties of the grid (cell areas and volumes, interpolation factors, normal 
distances of near-boundary cell nodes from boundary). All variables including grid node 
coordinates are converted to one-dimensional arrays. These are formed by scanning the grid in J- 
direction (figure A.2) for 1=1, and then repeating for all I's. The position of any node in one- 
dimensional array is therefore defined as; 

U = (I,J) = (I-l) * NJ + J 

The actual number of grid nodes is one row and one column less than for all cell nodes. For I = 
NI and J = NJ fictitious grid nodes are introduced which have the same coordinates as actual nodes 
on NI-1 in I-direction and NJ-1 in J-direction. 


- 24 - 



The subroutine then calculates interpolation factors which are associated with cell nodes and are 
used in the main program to calculate values of dependent variables at locations other than cell 
nodes (cell centers). Definition of these are given in Figure A.3. Cell areas and volumes are 
calculated next followed by calculations of normal distances of near-boundary nodes from all four 
outer boundaries. 


CALCE (PHI. IPHI) 

This subroutine solves the linearized and discretized transport equations for the turbulent energy k 
and the energy dissipation rate e. The two dummy parameters in the calling statement, PHI and 
IPHI, represent arrays containing dependent variables for which the equation is to be solved, the 
subroutine sets up the convective and diffusive coefficients over the entire field. Then it calculates 
the source terms for either k or £ transport equations. A call is made to entry MODPHI in order to 
modify these sources and boundary coefficients to suit the particular problem. Moreover, a check 
is made if the two-layer model is selected then the energy dissipation is set algebraically in the 
sublayer region. 

The discretized equations have the form 


A p 0 p = X A i + S 4> 

i=EWNS 

where the coefficients A,- ( i=E,W,N,S see figure A.3 ) contain both the convective and diffusive 
fluxes, these equations are assembled and solved by calling subroutine SOLSIP which is based on 
Stone’s Strongly Implicit Solver [19]. 


TWOLAY 


This subroutine is called if the two-layer or low-Reynolds number models are used. It calculates 
the different coefficients needed to describe the energy dissipation and eddy viscosity. In this 
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subroutine the normal distances used in the definition of the turbulent Reynolds number R y at 
comer regions are calculated based on the normal distance nearest to the solid boundary. 

SOLSIP 

This subroutine solves the system of linear algebraic equations for k and £ using Stone's Implicit 
Procedure [19]. The array RES (IJ) is used to store residuals. The sum of absolute residuals 
"RESORP" calculated in the first pass through this part of the routine is used as a measure of 
convergence of the solution process as a whole and this value is stored in RESOR (IPHI). This 
variable RESOR (IPHI) is passed to the main solver and if desired can be normalized and 
compared with the maximum error allowed there. If necessary, inner iterations counter L and the 
sum of absolute residuals RESORP are printed out to monitor the rate of convergence of k and £ 
solution. If the ratio RSM is greater than the maximum allowed for the variable in question, SOR 
(IPHI), and the number of inner iterations is smaller than a prescribed maximum, NSWP (IPHI), 
then the routine repeats the sequence of calculating the residuals, increment vectors and updating 
the dependent variable. 


USERM 

This subroutine has different ENTRY points or sections where variables are updated and boundary 
conditions are set. 

Section MODVIS 

This section calculates effective viscosity (Eq. 8). It is called after calculating k and £ . At locations 
where £ is close to zero (i.e., < 10 ' 30 ) viscosity is set to zero. A provision is made for under 
relaxing changes in effective viscosity which may help to stabilize oscillations and improve 
convergence rate. 


Section MODPHI 

This section is called from CALCE subroutine and sets the boundary conditions for k and £ 
depending on which variable being called (IDIR = 1 for k and IDIR = 2 for £). For the k -equation, 
the south boundary is checked first if it is one of four options: 
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(1) An inflow boundary ITBS(I) = 1, where the source term is set to accept the inlet values at 
J = 1 (south boundary) 

(2) Outflow boundary ITBS(I) = 2, where zero gradient in y or Ti-direction is employed. 

(3) Symmetry boundary, TBS(I) = 3, where gradients normal to symmetry plane are zero. 

(4) Wall boundary, ITBS(I) = 4, where the production term GENTS(I) calculated form 
subroutine WALLFN in program MODIFY is added to the rest of the source term 
SU(IJ). 

Boundary conditions for the £ -equation are similar to those of k except at the wall where they are 
set to appropriate values for each near wall treatment. 


A.4 Program MODIFY 

This program is compiled separately and is called from the u and v solver routines. It basically 
updates the flux source term of the discretized momentum equation due to wall shear stresses. If 
the u-momentum equation for example is discretized in the form 


A p u p 



* 

+ S u 


where P, E, W, N, S are cell nodes as shown in Figure A.3, and A p and A/'s contain convective 
and diffusive coefficients. S* is the source term containing pressure gradients and cross-derivative 

diffusion terms and convective terms for second-order upwinding scheme. This source term is 
usually linearized as 5* = S u - B p u p . The term B p is usually moved to the left hand side of the 

equation and modifies the diagonal coefficient A p = A* + B p , and the equation can be written as 


A P Un — y. At Ui + Su 

y i=EWNS 
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Then S u and B p are passed to subroutine MODIFY where they are modified if a wall is present 
(e.g., ITBS(I) = 4 for south boundary). 


For an iterative flow solver using the finite-volume methodology. A typical interface and call to the 
k-e module from the main flow solver can be represented by a flow chart as shown in figure A.4. 
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Figure A.l 2D/axisymmetric k-e module deck 
flow chart 
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In this section a description of the multi-time-scale k-e turbulence model that is coded as a self 
contained computer program to compute turbulent flow quantities in two-dimensional or 
axisymmetric geometry is given. Detailed description of the module structure, variables used and 
how to interface the module with CFD flow solvers are given in Appendix B. The module has been 
tested as a separate self-contained unit using the REACT code [1] and was independently tested at 
the University of Alabama at Huntsville (UAH) using own code (MAST). 


3.1 Introduction 

Turbulent flows comprise fluctuating motions with a spectrum of sizes and time scales and 
different turbulent interactions are associated with different parts of the spectrum. In the single- 
time-scale turbulence models such as the k-c turbulence model it is assumed that a single time scale 
(proportional to k/e) can be used to describe the turbulent flow. In many complex flows turbulence 
is generally in spectral inequilibrium and a single time scale description is a simplification. 

Figure 1, shows a sketch of a typical energy spectrum in a turbulent flow at high Reynolds number 
in a simplified split spectrum method. Two regions can be identified, the production range (at wave 
number K< k } ) where the kinetic energy (k p ) leaves this region at a raie (e p ) and a high wave 
number or dissipation region (k>Kj) with kinetic energy (k t ) and energy dissipation rate (£, ). 
Hanjalic et al. [2] developed a simple multiple-time-scale turbulence model based on a rational 
extension of the single scale equation ideas. In their model, a fixed ratio of the turbulent kinetic 
energy of large eddies ( k p ) to that of the fine scale eddies (k t ) is used to partition the spectrum. 
Kim and Chen [3] improved on the simplified split spectrum by dynamically determining the 
location of the partition (i.e kp/k t ) as part of the solution and is dependent on the turbulence 
intensity, production rate, energy transfer and dissipation rate. The variable partitioning method 
causes the effective eddy viscosity to decrease when production is high and to increase when 
production vanishes -a behavior consistent with experimental observations. 


3.2 Theory and Model Equations 

The multi-time-scale turbulence module is based on the variable partitioning of the turbulent energy 
spectrum proposed by Kim and Chen [3]. In this model the turbulent kinetic energy spectrum is 
divided into two sets of wave number regions giving two evolution equations for each region. 
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These equations represent the kinetic energy (k p ) and the energy transfer rate (£ p ) in the 
production range of the spectrum and the kinetic energy ( k t ) and the energy dissipation rate (e t ) in 
the dissipation range of the spectrum. This model allows the partition to move toward the high 
wave number region when production is high and toward the low wave number region when 
production vanishes. 

The equations which describe the multi-time-scale turbulence model used are given below. The 
turbulent kinetic energy and the energy transfer rate equations for the energy containing large 
eddies are given as; 



where G is the turbulence production rate, given as 



where fi is the viscosity 

Ht is the turbulent viscosity 

k p is the turbulent kinetic energy in the production range 

e p is the energy transfer rate 

ajcp and <Jsp are constants 

C p j , C p 2 and C p 3 are turbulent model constants 

The turbulent kinetic energy and the dissipation rate equations for the high wave number small 
scale eddies region are given as; 
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(4) 




t + pC « 



where k, is the turbulent kinetic energy in the dissipation range 
e, is the energy dissipation rate 
Ok, and o £t are constants 
C t j, C t 2 and C t 3 are turbulent model constants 

2 

7 £ 

The terms - C P j j~ and p C tl if represent variable energy transfer functions. The first term 

p k P kt 

increases the energy transfer rate when production is high and the second term increases the 
dissipation rate when the energy transfer rate is high. The turbulent viscosity is given as 

Ht - p Cpf - p Cfi 

£p £t 

where k= k p + k, is the total turbulent kinetic energy and Ctf is a constant. 

The model constants used are similar to those used by Kim and Chen [3] 

oic p = 0. 75, o £p - 1.15, Okt = 0.75, O et = 1.15 
Cpi = 0.21, Cp 2 = 1-24, Cp3 = 1.84, C tl = 0.29 
C t2 = 1.28, C t 3 = 1.66 and C pf = 0.09 


For turbulent flow analysis, equations (l)-(4) are solved by the module that is interfaced with a 
Reynolds averaged flow solver to compute the turbulent flow field. For an incompressible, steady 
and axisymmetric turbulent flow, a generalized equation that expresses the transport of turbulent 
flow can be written as; 
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where 0 is the dependent variable, which stands for 0 = u , v , w for the axial, radial and 
tangential velocities respectively, p is the fluid density, r& x and r®f are exchange coefficients in 
x and r -directions, respectively, and S® is the source term for the variable 0. 

The source terms for the dependent variable are: 

• Axial direction, 0= u, r<p x = 2p e , = p e and 

c dP 1 d ( dv 

Su = ~dx + 7 dF ( ^ e dx > 


where p e is the eddy viscosity and P is the pressure 

• Radial direction, 0= v, r<p x = p e , r<p r = 2p e and 

_ d ( du ) - v pw 2 dP 

Sv ~ ~ dx ' dr' 2 ^ e r 2 + r dr 


Tangential direction, 0- w, r<p x = p e , T® r - p e and 

pvw w d 


>W 


^ dr 
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Equations ( 1 )-(4) can also be written in a similar form as equation (5) where 0 stands for; 

• Turbulent kinetic energy in the production range of the energy spectrum 

Ut 

0= k p , r® = p + JZ -= /dv and 


^kn 
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Energy transfer rate in the production range of the energy spectrum 
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• Turbulent kinetic energy in the dissipation range of the energy spectrum 
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u, 

<P= k t , r<p Y - fl+ — = r<Ev and 
Okt 

Skt= P £ p - p £ t 

• Energy dissipation rate in the dissipation range of the energy spectrum 
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Near a wall, the wall function boundary conditions used axe similar to that of Kim and Chen [3]. 

A two layer model for the multi-time-scale k-£ turbulence model similar to that of Chen and Patel 
[4] for the single-time-scale k-£ turbulence model is included in the present release. 


3.3 Model Evaluation 

The multi-time-scale k-£ module was evaluated by comparisons with experimental studies. One of 
the test problems considered was the backward facing step of Driver and Seegmiller [5] where the 
multi scale k-£ model predicted a recirculation length of 6. 14 step heights (H) downstream of the 
step which is closer to the experimental value (6. 10 H) than the standard k-e model (5.35H). 

The majority of the tests were conducted using Roback and Johnson's experimental data [6] for 
swirling confined double concentric jets. Preliminary analysis indicated some sensitivity to the ratio 
kp/kt at the inlet boundary, however, a value of 3 was found reasonable in the present analysis. 
Figures 2 and 3 show the streamline patterns for wall functions and two-layer near wall treatments 
respectively. The upper (a) and lower (b) parts correspond to the single-scale k-e and the multi- 
scale k-e models respectively. It can be seen from these contours that there are two recirculation 
zones in the chamber, one is near the expansion comer and another located in the central region and 
accurate predictions of this central region is very important in combusting swirling flows. Figures 
4a and 4b, show the axial velocity along the centerline. In terms of strength and size of the central 
recirculation zone, the multi-scale k-e model yields better agreement than the single-scale k-e 
model. In the central recirculation region the k-e model tends to connect the energy transfer rate to 
the local mean strain rate too strongly while the multi-scale model suppresses this tendency. 
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Figures 5 and 6 show the radial profiles of the mean axial velocity at different axial locations 
downstream of the inlet using the wall function and the two-layer near-wall treatments respectively. 
Similarly, figures 7 and 8 show the corresponding profiles for the tangential velocity, and figures 9 
and 10 show the radial profiles of the axial normal turbulent intensity (mm ) m using both the wall 
function and the two-layer near-wall treatments. In general, the numerical results indicate that the 
multi-scale model gives better agreement than the standard k-£ model. 
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Figure 1. 



Description and nomenclature of the multiple-time-scale turbulence model; 


K=Ki 

k p = J E dK , k t = j E dK 

K=0 K=Kl 


ki = Partition wave number, E = Energy spectral density 
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Figure 2. Stream-function contours of confined swirling 

jet flow using the wall function near wall treatment 




(a) k-e model 



(b) M-S model 


Figure 3. Stream-function contours of confined swirling 
jet flow using the two-layer near wall treatment 
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Figure 4. Axial mean velocity along the centerline 
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(b) two-layer model 
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Figure 5. Radial profiles of mean axial velocity using 
the wall function near wall treatment 
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Figure 7 Radial profiles of mean tangential velocity 
using the wall function near wall treatment 
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Figure 9. Radial profiles of turbulent intensity (uu) 1 ' 2 
using the wall-function near wall treatment 
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Figure 10. Radial profiles of turbulent intensity 

( uu) 1 ' 2 using the two layer near-wall 
treatment 
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APPENDIX B 


Multi-Scale k-e Module Deck 


B.l Introduction 

This user's manual describes the multi-scale k-£ module deck. The module is a self contained 
FORTRAN source code to compute turbulent kinetic energy, energy dissipation and turbulent eddy 
viscosity using the multi-time-scale k-£ turbulence model. It uses as input the mean flow 
properties as computed by conventional CFD techniques. The module is constructed to be self- 
contained, stand alone and compatible with a number of CFD solvers. A discussion of the multi- 
time-scale k -£ module structure is given next together with flow charts to show how to interface 
the module with a number of flow solvers. A list of variable names used is also given. 


B.2 Program KEMOD 

This is basically the solver for the k and e - transport equations in both the production and the 
dissipation regions of the energy spectrum. It reads through its argument list different variables 
from the calling flow solver. These variables are described below where, each variable name ends 
with either an (I) for Integer variable, (R) for Real variable or (L) for Logical variable. 

The flow chart of the program is shown in Figure B.l. It shows the main operations performed by 
the code. 

List of Argument Variable Names 

XR Grid node locations in the x or ^-direction, dimensioned to XR (NX, NY) (input 

from flow solver) 

YR Grid node locations in the y or ri-direction, dimensioned to YR (NX, NY) (input 

from flow solver) 
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UR Axial or x-direction velocity (u), dimensioned as UR (NX, NY) (input from 

flow solver 

VR Radial or y-direction velocity (v), also dimensional as VR (NX, NY) (input 

from flow solver) 

WR Azimuthal velocity (w), dimensional WR (NX,NY) (input from flow solver) 

TER Large scale turbulence kinetic energy k p , dimensioned TER (NX,NY) 

(calculated in KEMOD and returned to flow solver) 

EDR Large scale turbulent energy dissipation rate £ p , dimensioned EDR (NX, NY) 

(calculated in KEMOD and returned to flow solver) 

TETR Small scale turbulence kinetic energy k t , dimensioned TETR (NX,NY) 

(calculated in KEMOD and returned to flow solver) 

EDTR Small scale turbulent energy dissipation rate e f , dimensioned EDTR (NX,NY) 

(calculated in KEMOD and returned to flow solver) 

DENR Fluid density, dimensioned DENR (NX,NY) 

URFKER Under-relaxation factors dimensioned as URFKER(4) and specified as follows: 
URFKER(l) for large scale turbulent energy equation 
URFKER(2) for small scale turbulent energy equation 
URFKER(3) for large scale turbulent energy dissipation equation 
URFKER(4) for small scale turbulent energy dissipation equation 

PRTKER Prandtl/Schmidt numbers dimensioned as PRTKER(4) and specified as 
follows: 

PRTKER(l) for large scale turbulent energy equation 
PRTKER(2) for small scale turbulent energy equation 
PRTKER(3) for large scale turbulent energy dissipation equation 
PRTKER(4) for small scale turbulent energy dissipation equation 

GR = 1 .0 if second order upwinding is desired 
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= 0.0 if first order upwinding is used (input from flow solver). 

FIR Mass flux variable at cell faces in x- or ^-direction, dimensioned FIR (NX,NY) 

(input from flow solver) 

F2R Mass flux variable at cell faces in y or redirection, dimensioned F2R (NX, NY) 

(input from flow solver) 

UERI Iteration number (input from flow solver), this number must be equal to 1 for a 

restart case 

VISCOSR Dynamic viscosity (input from flow solver) 

VISR Eddy viscosity, dimensioned VISR (NX,NY) (calculated in KEMOD and 

returned to main solver) 

URFVISR Under-relaxation factor for total viscosity calculation 

AKSIL Logical variable for axisymmetric geometry (AKSIL=-TRUE-) or plain 

geometry (AKSEL=- FALSE-) (input from flow solver) 

C1R Turbulence model constant, C\ (input from flow solver) 

C2R Turbulence model constant, C 2 (input from flow solver) 

CMUR Turbulence model constant, (input from flow solver) 

NIMI Number of cell nodes in the I- or ^-coordinate lines, (input from flow solver) 

NJMI Number of cell nodes in the J- or r| -coordinate lines, (input from flow solver) 

JTBEI Boundary condition flag along east boundary must have one for each boundary 

node set to: 1-inlet, 2-outlet, 3-symmetry and 4-wall e.g., for an outlet 

boundary condition on the east boundary set JTBEI to NJ*2, and similarly for 
other boundaries, dimensioned JTBEI (NY) (input from flow solver) 
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JTBW1 Boundary condition flag along west boundary, dimensioned JTBWI (NY) 

(input from flow solver) 

ITBNI Boundary condition flag along north boundary, dimensioned ITBNI (NX) 

(input from flow solver) 

ITBSI Boundary condition flag along south boundary, dimensioned ITBSI (NY) 

(input from flow solver) 

Program KEMOD is interfaced with the main flow solver by a call to KEMOD with its arguments. 
For iterative flow solvers KEMOD is called within the iteration sequence after the solution of the 
momentum equations where the mean velocities are passed to KEMOD. There are different flow 
solvers utilizing different schemes from staggered to nonstaggered grid arrangement and for 
nonorthogonal coordinate system there are at least three alternatives to the choice of the velocity 
components 

i. Cartesian velocity components 

ii. Contravariant velocity components 

iii. Covariant velocity components 

The Cartesian velocity components are the most widely used and have the advantage of simple 
formulation of the governing equations. Whatever the arrangement used, mass fluxes at cell faces 
are required and passed to KEMOD as FIR and F2R in both directions. The location of other 
variables such as k and £ are at the cell center or cell nodes. 

The module starts by reassigning variable names passed to it from flow solver to names that are 
shared with the different subroutines of the module in an include file "mske.h". The user must set 
the values for NX and NY in mske.h greater than or equal to the maximum grid dimensions. Then 
a check is made if it is the first iteration in which case the grid file "GRIDG" is called -after passing 
the grid node locations XR & YR in KEMOD- in order to calculate grid related quantities which 
will be explained later. The need to call GRIDG can be waived if all the grid data are passed to the 
module. That is all the information about the grid such as interpolation factors FX and FY, cell 
areas (ARE) and volumes (VOL) and normal distances of first grid point from grid boundaries 
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(DNS from south boundary, DNN - from north boundary, DNW - from west boundary and DNE - 
from east boundary). 

After this, two calls to subroutine CALCKE are made to calculate the large and small scale 
turbulent kinetic energies with the identifier IPHI=1 and 2 respectively). The large and small scale 
energy dissipation equations are solved next by calling subroutine CALCKE again with the 
identifiers IPHI=3 and 4 respectively. The effective viscosity is calculated next. At locations 
where £ is close to zero (i.e., < 10 30 ) viscosity is set to zero. A provision is made for under 
relaxing changes in effective viscosity which may help to stabilize oscillations and improve 
convergence rate. 


B.3 Subroutines 
GRIDG 

Before calling this subroutine, the coordinates of all grid nodes, defined in reference to a fixed 
Cartesian coordinate frame are read. Figure B.2 shows the position of cell and grid nodes. 

This subroutine is called only once to calculate coordinates of grid nodes (intersection of grid lines) 
and geometrical properties of the grid (cell areas and volumes, interpolation factors, normal 
distances of near-boundary cell nodes from boundary). All variables including grid node 
coordinates are converted to one-dimensional arrays. These are formed by scanning the grid in J- 
direction (figure B.2) for 1=1, and then repeating for all I's. The position of any node in one- 
dimensional array is therefore defined as; 

1J = (1,J) = (1-1) * NJ + J 

the actual number of grid nodes is one row and one column less than for all cell nodes. For I = NI 
and J = NJ fictitious grid nodes are introduced which have the same coordinates as actual nodes on 
NI-1 in I-direction and NJ-1 in J-direction. 

The subroutine then calculates interpolation factors which are associated with cell nodes and are 
used in the main program to calculate values of dependent variables at locations other than cell 
nodes (cell centers). Definition of these are given in Figure B.3. Cell areas and volumes are 
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calculated next followed by calculations of normal distances of near-boundary nodes from all four 
outer boundaries. 


CALCKE fPHI. IPHI) 

This subroutine solves the linearized and discretized transport equations for the turbulent energies 
(k p and k t ) and the energy dissipation (£ p and e t ). The two dummy parameters in the calling 
statement, PHI and IPHI, represent arrays containing dependent variables for which the equation is 
to be solved, the subroutine sets up the convective and diffusive coefficients over the entire field. 
Then it calculates the source terms for either k or e transport equations. A call is made to entry 
MODMSKE in order to modify these sources and boundary coefficients to suit the particular 
problem. 

The discretized equations have the form 


A p <P p = £a< + S 4> 

i=EWNS 

where the coefficients A, ( i=E,W,N,S see figure B.3) contain both the convective and diffusive 
fluxes, these equations are assembled and solved by calling subroutine SOLSIP which is based on 
Stone’s Strongly Implicit Solver [7]. 


SOLSIP 

This subroutine solves the system of linear algebraic equations for k and £ using Stone's Implicit 
Procedure [7]. The array RES (IJ) is used to store residuals. The sum of absolute residuals 
"RESORP" calculated in the first pass through this part of the routine is used as a measure of 
convergence of the solution process as a whole and this value is stored in RES OR (IPHI). This 
variable RESOR (IPHI) is passed to the main solver and if desired can be normalized and 
compared with the maximum error allowed there. If necessary inner iterations counter L and the 
sum of absolute residuals RESORP are printed out to monitor the rate of convergence of k and e 
solution. If the ratio RSM is greater than the maximum allowed for the variable in question, SOR 
(IPHI), and the number of inner iterations is smaller than a prescribed maximum, NSWP (IPHI), 
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then the routine repeats the sequence of calculating the residuals, increment vectors and updating 
the dependent variable. 


MODMSKE 

This subroutine is called from CALCKE subroutine and sets the boundary conditions for kp, kj 
and £p, £t depending on which variable being called (EDIR = 1, 2, 3, and 4 for kp, k t , £ p , and £ t 
respectively). Consider the south boundary for example, if it is one of four options: 

(1) An inflow boundary ITBS(I) = 1, where the source term is set to accept the inlet values at 
J = 1 (south boundary) 

(2) Outflow boundary ITBS(I) = 2, where zero gradient in y or T|-direction is employed. 

(3) Symmetry boundary, TBS(I) = 3, where gradients normal to symmetry plane are zero. 

(4) Wall boundary, ITBS(I) = 4, where the production term GENTS(I) calculated form 
subroutine WALLFN in program MODIFY is added to the rest of the source term 
SU(IJ). 


B.4 Program MODIFY 

This subroutine is called from the u and v solver routines. It basically updates the flux source term 
of the discretized momentum equation due to wall shear stresses. If the u-momentum equation for 
example is discretized in the form 


A p u p 


I Ai 

i=EWNS 


Ui 


* 

+ S u 


where P, E, W, N, S are cell nodes as shown in Figure B.3, and A* and A, 's contain convective 
and diffusive coefficients. 5* is the source term containing pressure gradients and cross-derivative 
diffusion terms and convective terms for second-order upwinding scheme. This source term is 
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usually linearized as S u — S u - BpUp The term Bp is usually moved to the left hand side of the 
equation and modifies the diagonal coefficient A p = A* + B p , and the equation can be written as 

Ap u„ = X + Su 

y i=EWNS 

Then S u and B p are passed to subroutine MODIFY where they are modified if a wall is present 
(e.g., ITBS(I) = 4 for south boundary). 

List of Argument Variable Names 

CMU Turbulence model constant, (input from flow solver) 

VISCOS Dynamic viscosity (input from flow solver) 

XX Grid node locations in the x or ^-direction, dimensioned to XX (NX*NY) 

(input from flow solver) 

YY Grid node locations in the y or r\ -direction, dimensioned to YY (NX* NY) 

(input from flow solver) 

R Grid node radius equal to 1 for non-axisymmetric and YY for axisymmetric, 

dimensioned to R (NX*NY) (input from flow solver) 

DNS Normal distance to south, dimensioned to DNS (NX*NY) (input from flow 

solver) 

DNN Normal distance to north, dimensioned to DNN (NX*NY) (input from flow 

solver) 

DNW Normal distance to west, dimensioned to DNW (NX*NY) (input from flow 

solver) 

DNE Normal distance to east, dimensioned to DNE (NX*NY) (input from flow 

solver) 
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U Axial or x-direction velocity (u), dimensioned as UR (NX*NY) (input from 

flow solver 

V Radial or y-direction velocity (v), also dimensional as VR (NX*NY) (input 

from flow solver) 

W Azimuthal velocity (w), dimensional WR (NX*NY) (input from flow solver) 

DEN Fluid density, dimensional DEN (NX* NY) (input from flow solver) 

TE Large scale turbulence kinetic energy k p , dimensioned TE (NX*NY) 

(calculated in KEMOD and returned to flow solver) 

TET Small scale turbulence kinetic energy k t , dimensioned TET (NX*NY) 

(calculated in KEMOD and returned to flow solver) 

SU Variable source term, dimensioned SU (NX*NY) 

BP Constant source term, dimensioned BP (NX*NY) 

AE Cell area, dimensioned to AE (NX*NY) (input from flow solver) 

AW Cell area, dimensioned to AW (NX*NY) (input from flow solver) 

AN Cell area, dimensioned to AN (NX*NY) (input from flow solver) 

AS Cell area, dimensioned to AS (NX*NY) (input from flow solver) 

SUVS,SPVS,SUWS,SPWS 

Source terms at south boundary due to wall shear stress, all dimensioned to 
SMS (NX*NY) (returned to flow solver) 

SUVN,SPVN,SUWN,SPWN 

Source terms at north boundary due to wall shear stress, all dimensioned to 
S##N (NX*NY) (returned to flow solver) 
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SUVW,SPVW,SUWW,SPWW 

Source terms at west boundary due to wall shear stress, all dimensioned to 
S##W (NX*NY) (returned to flow solver) 

SUVE,SPVE,SUWE,SPWE 

Source terms at east boundary due to wall shear stress, all dimensioned to S##E 
(NX*NY) (returned to flow solver) 

GENTS ,GENTN ,GENTW , GENTEE 

Generation terms at south, north, west, and east boundaries respectively due to 
moving walls, with GENTS(NX), GENTN(NX), GENTW(NY), and 
GENTEE(NY) (returned to flow solver) 

NX Maximum number of cell nodes in the I- or ^-coordinate lines, (input from flow 

solver) 

NY Maximum number of cell nodes in the J- or -coordinate lines, (input from flow 

solver) 

NXNY NX*NY 

NIM Number of cell nodes in the I- or ^-coordinate lines, (input from flow solver) 

NJM Number of cell nodes in the J- or T|-coordinate lines, (input from flow solver) 

ITBS Boundary condition flag along south boundary must have one for each 

boundary node set to: 1 -inlet, 2-outlet, 3-symmetry and 4-wall e.g., for an 
outlet boundary condition on the east boundary set ITBS to NI*2, and similarly 
for other boundaries, dimensioned ITBS (NX) (input from flow solver) 

ITBN Boundary condition flag along north boundary, dimensioned ITBN (NX) (input 

from flow solver) 

JTBW Boundary condition flag along west boundary, dimensioned ITBNI (NY) (input 

from flow solver) 
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JTBE 


Boundary condition flag along east boundary, dimensioned JTBE (NY) (input 
from flow solver) 


For an iterative flow solver using the finite-volume methodology. A typical interface and call to 
KEMOD from the main flow solver can be represented by a flow chart as shown in figure B.4. 
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Figure B.l Multi-scale k-e module deck flow chart 
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Figure B.4 Typical main flow solver flow chart with 
calls to the multi-scale k-e module 
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4.1 Introduction 


In this section a description is given of the two-dimensional/ Axisymmetric Algebraic Stress 
turbulence Model (ASM) based on the work of Rodi [1]. The model is coded as a self contained 
computer program to compute turbulent flow quantities when interfaced with a CFD solver. 
Detailed description of the module structure, variables used and how to interface the module with 
CFD flow solvers are given in Appendix C. 

The module uses as input the mean flow properties, as computed by conventional CFD solvers, 
and calculates the Reynolds stresses, turbulent kinetic energy and the energy dissipation. It is 
structured to be self-contained and compatible with many CFD codes. It has been tested as a 
separate unit at Rocketdyne using the finite-volume REACT code [2], The module has also been 
tested independently at the University of Alabama at Huntsville (UAH) using own code MAST. 

The module computes turbulent flow quantities in two-dimensional planar or axisymmetric 
geometry with or without swirl. The standard wall functions and the two-layer model of Chen and 
Patel [3] are used for the near wall treatment. 


4.2 Theory and Model Equations 

The Algebraic Stress (ASM) module is based on the work of Rodi [1]. The idea is to simplify or 
truncate the Reynolds stress equation by approximating the convective and diffusive transport of 

the Reynolds stresses uiu] in terms of the corresponding transport of turbulent energy. This 
allows the transport equation for the stresses to be expressed as a set of algebraic formulae 
containing the turbulence energy and its rate of dissipation as unknowns in the form: 

k 2 e 

UiUj = JpT) [ Pij ' 3 5ij £ + 0ij} 

where Pij = Production and P =2 Pkk and 

&ij = pressure-strain redistribution 

<P i'j = 4> jj t i + 4 >jj'2 + 4>ij,1w + 4>jj'2w 


- 49 - 



Rotta’s linear retum-to-isotropy concept for the non-linear part of 

e 2 

®ij,1 = -Ci £ (uiuj - 3 k Sjj) 

is used and the "isotropization of production" concept for the linear "rapid" part of 

®ij,2 = -C 2 ( Pi j - | P Sij) 

is used. Gibson and Launder [4] concept for the wall reflection terms is used as 

e c 3 3_ 

<Pij,lw = C] w pj( UkU m nkn m Oij - ^ “Wi nknj - 3 u k u j n kni ) J 

3 3 

®ij,2w = C2w ( ®km,2 nkn m 8ij - J &ik,2 nkrij - 2 ®jk,2 n^i )f 


where (np is the wall-normal unit vector in the i -direction. The wall-distance function (f) 

k 3/2 

represents the ratio of the turbulence length scale ( L e = ) and the wall distance and is given 


as 


cf 5 *' 5 1 

f , -JH j 

} ( K e ’An 

with An being the wall-normal distance. 


The set of algebraic stress equations can be arranged in the form 

Ay «2 + Bij v2 + Cij + Dij MV + Eij vw + Fy uw= Gij 

where Ajj, Bjj, Cij, Djj, Eij, Fjj, and Gij are functions of the mean and turbulent flow 
variables. 


The above equation can be solved iteratively in the main flow solver. However, the algebraic 
system of equations is stiff and convergence difficulties are encountered when solved iteratively. 
Therefore, the set of equations was cast in the general matrix form A T=B, where 
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3e _ 
2Ax- + 2 

dU 

dx 

dV 
' dy 

V 

r 

*¥y 

dV 

dx 

dW W 

' dy + r 

dW 

dx 

d_U 

dx 


3 ^ + 2 ^ 
2Xk dy 

V 

r 


dJJ 

dy 

,dW n W, 
1 dy r ' 

dW 

1 ' dx 

dU 


dV 

3e „ V 

dU 

dV. 

dW W 

0 dW 

dx 


' dy 

2Xk + ^ r 

"* dy + dx 

2 + 
dy r 


dV 

dx 


dU 

dy 

0 

_£_ dJJ dV 
A k + dx + dy 

0 

W 

r 



dW 

W 

dW 


dV V 

dV 

0 


dy 

r 

dx 


Xk + dy + r 

dx 

dW 



0 

dW 


dU 

_£_ dLJ 

dx 


0 

dy 


dy 

A k + dx + 


7 = [puu,pvv,pww,puv,pvw,puw ] 

p£ 3 . ^ . 

B = J + 2(UCi) ' 011 ’ 1w + 0 l^2w) 

p£ 3 , 

X + 2(1 -C 2 ) ( 022 ’ 1w + <p22 ’ 2w ' 
p£ 3 ^ , 

X + 2(1 -C 2 ) (® 33 ’ 1w + ® 33 > 2w ' 

(&12,1w + &12,2w) 

(1-C 2 ) (® 23 ’ 1w + ® 23 ’ 2w ) 

(I-C2) (® 13 ’ 1w + ® 13 > 2w ) 

1-C2 

where X = □ 

Ci-1+ ~ 

P £ 
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The matrix was inverted at each iteration step to obtain a converged solution. The wall function and 
a two-layer model were built in the module to model the near-wall region. 


4.3 Module Evaluation 

The ASM module was evaluated by comparison with experimental data of Driver and Seegmiller 
[5] for the backward facing step and the data of Roback and Johnson [6]. The effect of the wall 
reflection term is also studied with both wall function and two-layer near wall models. Figures la 
and lb show the stream-function contours for a backward facing step flow using the wall function 
and the two-layer near wall models with reattachement length of 5.59H and 5.83H respectively (H 
is the step height). Figures 2a and 2b show the stream-function contours for the Roback & 
Johnson confined swirling jet flow using the wall function near wall model, where (a) includes the 
wall-reflection term in the pressure-strain redistribution term and (b) without the wall reflection 
term. Figure 3 shows a comparison of the axial velocity along the centerline with and without wall 
reflection term. The comparisons of the predicted mean axial velocity, mean tangential velocity, 

turbulent intensities , vv , ww^ , and the Reynolds stress uw using the ASM model as 
compared with the single and multi-scale k-e models are presented in figures 4 to 9 respectively. 
The figures in general show that the ASM model used here when combined with the wall function 
near wall treatment predicts better comparisons without using the wall reflection terms. This may 
be explained by the fact that the wall reflection terms -whose purpose is to damp normal turbulent 
intensity normal to the wall as the wall is reached- are not effective when using wall functions near 
the wall. Similar conclusions were also obtained by the UAH group when testing the ASM module 
using their code (MAST). Also, in the ASM model, a set of algebraic equations for the Reynolds 
stresses are solved and there is no boundary conditions are needed for the stresses. This is not the 
same in the full Reynolds stress model (RSM) where a set of nonlinear differential equations are 
solved and boundary conditions for the stresses are required. More on this will be discussed in 
detail in the next RSM module. Also, more details will be given on the tensorial incorporation of 
the wall reflection terms since they are tied to the orientation of the wall through the unit normal 
vectors. 
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(a) 



(b) 


Figure 1. Stream-function contours of backward facing step 
flow using the ASM with (a) wall function and 
(b) two-layer near wall treatment 
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Figure 2. Stream-function contours of confined swirling 
jet flow 
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Figure 3. Decay of axial velocity along the centerline 
in confined swirling jet flow 
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Figure 5. Radial profiles of mean tangential velocity 
in confined swirling jet flow 
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Figure 7. Radial profiles of turbulent intensity (vv) i/2 
in confined swirling jet flow 
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Figure 8. Radial profiles of turbulent intensity (ww) m 
in confined swirling jet flow 
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APPENDIX C 


2D/Axisymmetric Algebraic Stress Turbulence Module Deck 

This module is a FORTRAN source code to solve 2D/Axisymmetric turbulent flow quantities using 
the algebraic stress model when interfaced with a main flow solver. The module consists of the 
main routine ASMOD that calls a number of subroutines to perform different functions that will be 
explained below. 

3 . 1 Subroutine ASMOD 

This is basically the main routine that reads through its argument list different variables from the 
calling flow solver which are described below. 

List of Argument Variable Names 

X Grid node locations in the x or ^-direction, dimensioned to X(NX*NY) 

Y Grid node locations in the y or ri-direction, dimensioned to Y(NX*NY) 

FX Interpolation factor in the x or ^-direction. 

FY Interpolation factor in the y or ri-direction. 

ARE Cell areas 

VOL Cell volumes. 

R Radial distance in the axisymmetric geometry or 1. for planar geometry. 

DNS Normal distance of a cell from the south-boundary dimensioned to NX. 

DNN Normal distance of a cell from the north-boundary dimensioned to NX. 

DNE Normal distance of a cell from the east-boundary dimensioned to NY. 

DNW Normal distance of a cell from the west-boundary dimensioned to NY. 

U Axial or x-direction velocity, dimensioned to NX*NY. 

V Radial or y-direction velocity, dimensioned to NX*NY. 

W Tangential or azimuthal velocity, dimensioned to NX*NY. 

TE Turbulent kinetic energy, dimensioned to NX*NY. 

FT) Turbulent energy dissipation rate, dimensioned to NX*NY. 

DEN Density (assumed constant for incompressible flows). 

FI Mass flux at cell faces in the x or ^-direction, dimensioned to NX*NY. 

F2 Mass flux at cell faces in the y or r^-direction, dimensioned to NX*NY. 
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VISCOUS Laminar viscosity. 

VIS Eddy viscosity, dimensioned to NX*NY. 

RESOR Residual error for the k and £ -equations solver, dimensioned to 2. 

ITBS Boundary condition flag along the south boundary dimensioned to NX and 

must have one for each boundary node set to: 1 -inlet, 2-outlet, 3-symmetry 
and 4-wall, e.g., for a wall boundary condition along the south boundary 
set ITBS to NX*4. Similarly for the other boundaries. 

UBN Boundary condition flag along the north boundary, dimensioned to NX. 

JTBE Boundary condition flag along the east-boundary dimensioned to NY. 

JTBW Boundary condition flag along the west-boundary dimensioned to NY. 

ITER Iteration number. 

FMUU function used in the two-layer model. 

ICAL = 1 for swirl velocity calculations, 0 otherwise. 

AKSI = 1 for axisymmetric flow, 0 otherwise. 

RESTART = 1 if calculations are restarted from a previous run, 0 otherwise. 

ASMOD starts by reading the turbulent flow constants, under-relaxation factors and 

Prandtl/Schmidt numbers for the k and £ equations. These are; 

CD1, CD2 constants in the k and £ -equations and are usually set to 1.44 and 1.92 

respectively. 

CMU, ELOG, constants in the k and £ -equations and are usually set to 0.09, 9.8 and 0.42 

and CAPPA respectively. 

LAY2 set to true (T) for two-layer model and false (F) for wall functions. 

GKE is set to 1 for second-order upwinding of the convective terms in the k and 

£ -equations. 

ALFAKE is the iteration parameter used in the k and £ -equation solver. 

URFVIS is the underrelaxation factor of the viscosity near the wall. 

SORKE( 1 ) and are the degree of accuracy for the k and £ -equation solver respectively. 

SORKE (2) 

URFKE( 1 ) and are the underrelaxation factors for the k and £ -equations respectively. 
URFKE(2) 

PRTKE(l) and are ratio of Prandtl to Schmidt numbers used in the k and £ -equations in the 

PRTKE(2) two-layer model near the wall. 

C1,C2 are constants in the ASM model. 
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C1W and C2W 


are the two constants in the wall-reflection terms of the pressure-strain 
redistribution term. 

CK and CE constants in the diffusion term of the k and £ -equations. 

WREFON = 1 if the wall reflection terms of the pressure-strain term are to be 

included, o otherwise. 

All dimensions considered are one-dimensional. The position of any node is defined as IJ = (I,J) = 
(I-1)*NJ +J, where NI and NJ are the number of grid nodes in the X and T-directions 
respectively. It is assumed that grid related data such as cell areas, volumes and interpolation 
factors be passed to the module from an external grid generator. 


Subroutine WALREF 

This subroutine calculates the wall reflection terms in the pressure-strain redistribution term. It 
calculates the wall unit normal vectors and the normal distance away from the wall. This is needed 
to resolve the wall tangential and normal velocity components that are needed to obtain the near- 
wall values of the Reynolds stresses. 

Subroutine CALPIJ 

This subroutine calculates the production terms of the individual stress components. 

Subroutine CALUIUJ 

This subroutine calculates the individual stress component from its algebraic equation. It sets the 
coefficients of the algebraic stress equations which are solved implicitly at each iteration step by 
inverting a 6x6 matrix. 

Subroutine ACALCKE 

This subroutine solves the transport equations for the turbulent energy (IPHI=1) and energy 
dissipation.(IPHI=2). Daly and Harlow [7] gradient stress diffusion form is used in the module 
instead of the simplified isotropic diffusivity form. The subroutine calls MODPHI subroutine that 
sets the appropriate boundary conditions for k and £ . The set of algebraic difference equations are 
then solved using Stone's strongly implicit solver ASOLSIP. 
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Subroutine ATWOLAY 


This subroutine calculates the near wall turbulence using Chen and Patel’s [3] two layer model. 

Subroutine MODPIJ 

This subroutine modifies the production terms near the wall using the near wall region model. 

Subroutine MODPHI 

This subroutine calculates the near wall boundary conditions for the turbulence energy and the 
energy dissipation. 

Subroutine AMODVIS 

This subroutine modifies the eddy viscosity close to a wall using the near wall model chosen. 

Subroutine ASOLSIP 

This subroutine solves the system of linear algebraic equations for k and e using Stone's Implicit 
Procedure [8], 

Subroutine AMODIFY 

This subroutine is called from the momentum equations solver of the main routine. It updates the 
flux source terms of the discretized momentum equations due to wall shear stresses and due to the 
Reynolds stress gradients. The terms SUASM, SVASM and SWASM need to be added to the U, 
V and W-momentum equations of the main solver. They represent the difference form of the 
Reynold stress gradients in the momentum equations. 
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5.1 Introduction 


This report describes a self contained FORTRAN source code to compute turbulent quantities 
using Launder, Reece and Rodi’s [ 1 ] second order closure, Reynolds stress model The module 
deck is designed to interface with a number of flow solvers to analyse incompressible turbulent 
internal flows. Detailed description of the model used is given with a special emphasis on the 
coupling of the mean velocity and Reynolds stresses in the discretization procedure of the 
generalized coordinate system using a co-located finite volume method. The module was interfaced 
with the REACT flow solver and tested with benchmark flows including the backward-facing step. 
The module was also successfully interfaced with the MAST code at the University of Alabama at 
Huntsville (UAH) and independently tested. The Reynolds stress model implemented produced 
consistently more accurate simulations that the standard k-e model. 


5.2 Theory and Model Equations 


The flow is considered planar or axially symmetric, steady with constant fluid properties. Its mean 
field may be described by a two-dimensional time averaged equations of continuity and 
momentum, which can be written as; 


dpU i dprV . 

~d^ + T ~dT = 0 


( 1 ) 


d(prU<P) d(prV<P) d . d<P d , d<P 

^ r - + ~ = s 1 < r ^> * s < r ^> + rS ° 


(2) 


<P stands for any of the dependent variables, namely, U and V (axial and radial velocities 
respectively) and rW (radial distance r multiplied by the tangential velocity W ). p is the fluid 
density, p is the laminar viscosity. p is the source term for the variable <J> and is given by; 


dP dpu 2 1 dpruv 

- Axial direction, <P = U and Sy = - ^ ^ ^ 


- Radial direction, 0-V and Sy = 


dP pW 2 2pV ]_ drrv 2 druv rw 2 
dr + r r 2~r dr dx + r 


_ „ r . „ „ p drW ruw drvw „ — 

- Tangential direction, <P - rW and Sw = - 2 — - p - - p ^ - 2rvw. 
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where u, v and w are the fluctuating velocity components in the axial, radial and azimuthal 
directions respectively. 

Turbulence wall effects in the module are represented by Gibson and Launder [2] version of the 
high Reynolds number stress transport closure of Launder, Reece and Rodi [1]. The stress closure 

consists essentially of modeled transport equations for the stresses UjUj and for axisymmetric 
swirling flow it includes all the six stresses u 2 , v 2 , w 2 , uv, uw and vw. 

The set of differential equations governing the transport of Reynolds stresses (-u^j ) is obtained 
from Navier-Stokes equations by multiplying the equations for the fluctuating components («,• ) 
and (uj ) by (uj ) and («,• ) respectively, then summing these equations and time averaging the 
results. The resulting Reynolds stress transport equations are then solved using the mean flow 
equations to obtain the mean and turbulent flow quantities. 

The full transport equations for the Reynolds stresses can be written in a compact form using 
Cartesian tensor representation as; 


1 dprUjc UiU[ 
r dxk 


rdxk 


(rCkP ukui- 

e 


d UjUj 
dxi 


) = Pij + Dij + 




£ ij 


(3) 


Where Uk are the mean velocity components in xk -direction. The right hand side contains the 
production term Py given as 

„ , — dUj — dUi , 

Pij = -P(Wk-££+u jUk -^) (4) 

Pij does not require approximations since it is fully represented by turbulent stresses and mean 
flow gradients. 

The dissipation correlation £y arise from the fine-scale of the turbulent motion. At high Reynolds 
numbers these scales are many orders of magnitude smaller than the large energy containing eddies 
and turbulence energy cascades down along the eddy-size range with little linkage occurring at 
intermediate scales, to be ultimately dissipated by the smallest eddies which are unaware of the 
nature of the mean flow and the large scale turbulence. Therefore, the structure of these fine scale 
motions responsible for viscous dissipation is isotropic and the dissipation tensor reduces to 
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(5) 


2v $5trl cS » 


An additional equations for the dissipation £ is required. 


Dij represents the Reynolds stress diffusion which does not in general contribute greatly to the 
balance of transport of uiiij except in regions of low stress production by mean strain. This term 
include contributions of fluctuating pressure-velocity correlations ( put and puj ), triple correlations 

uiujuk and viscous diffusion v . Daly and Harlow [3] proposed a simple gradient diffusion 
hypothesis to model the stress diffusion term in the form 


Di J ~ Cs dx k 


k , — duiUj , 

p - e Ium W 


(6) 


with constant C s is taken to be 0.22. Lien and Leschziner [4] simplified the treatment of the 
diffusion term to allow an appropriate isotropic diffusivity in the form 


,. = J_ r }LJL 

tJ dx k 1 Gk dx k 


(UiUj)] 


(7) 


where <J k is a dimensionless constant. Harlow's proposal for the diffusion term is adopted in the 
present module since it is based on the fundamental conservation equations for the triple 
correlations, while Lien & Leschziner's form has a weaker basis in this respect. 

represents the redistribution of turbulence energy among the normal stresses through the 
interaction of pressure and strain fluctuations. Modeling the pressure-strain term is the most 
elaborate and involves the solution of the Poisson equation for pressure fluctuations p. The explicit 
appearance of the pressure in the correlation is eliminated by taking the divergence of the equation 
for the fluctuating velocity m, , thus obtaining a Poisson equation for p. Following a volume 
integration of the resulting equation subject to the assumption of local mean-flow homogeneity 
results in three contributions to the pressure-strain correlation <£y . One involving just fluctuating 
quantities & lJ} ] another arising from the presence of the mean rate of strain <Pij,2- and a third 
arising from the surface integral representing wall effects d>y>. Since the primary role of <f>y is to 
guide turbulence towards isotropy, Rotta [5] proposed for ; 

<Pijj = - 2 p C] ebij (8) 
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— 2 

where bij = ( ujuj - j <5y / k) /2k is the dimensionless anisotropy parameter. C/ is a constant and A: 

and e are turbulent kinetic energy and energy dissipation respectively. More elaborate models have 
been proposed such as Lumley [6] and Fu [7] using a nonlinear expression for <%;. The term 
<Pij 2 has been the subject of more extensive research. The traditional linear approach similar to 
Rotta's work simplifies this correlation to: 

®ij,2 = -C 2 (Pij- 1 SijP) (9) 

where P is the production of turbulent kinetic energy. Analogous to <Pijj the correlation, &ij : 2 
represents the isotropization of turbulence production tensor with C 2 as a constant. More elaborate 
models such as that of Speziale, Sarkar and Gatski [8] is based on dynamical systems approach 
and invariancy concepts. Nonlinear models for 4^2 based on the realizability constraints have 
been developed, e.g, Shih and Lumley [9] and Fu, Launder and Tselepidakis [10], The simplified 
correlations in equations (8) and (9) are used in the present module. 


The correlation 0 t j iW represents the wall damping effects that counteracts the tendency of dfy i 
and <Py,2 to isotropise the turbulent structure. Since close to a solid wall turbulence approach a 
state of intense anisotropy associated with a tendency towards a 2D turbulence. Following Shir 
[11] and Gibson and Launder [2], &ij,w is modeled as the combination of two separate terms; 

£ 3 3 i 

0[jJw — C] w p ^ [ UfcUm n k n m &ij ~ 2 u k u i n k n j ~ 2 u k u j n k n i 7 f ( J^) (10) 

3 3 l 

0ij_2w = C 2 w [ ^km.l n k n m &ij * 2 n k n j ~ 2 ^jk,2 n k n i 1 f ( J^) (11) 


i tf/2 

where l„ is the normal distance from the point in question to the wall and l ( = ) is the turbulent 

e 


length scale. The following relationship is used for the wall damping function 
/ = 


Cl 5 #" 1 


( 12 ) 


k e <ln> 

where </„> is the average distance of the point considered from the surrounding surfaces and n, is 
a wall-normal unit vector in the i -direction. The constants Cj w and C 2w have values of 0.5 and 
0.3 respectively. 


It will be of some value to list the full Reynolds stress equations for axisymmetric swirling flows. 
Although, the derivations have been carried out within the constraints of Cartesian coordinates, 
considerations will be given next to the forms applicable to any general curved coordinate system. 
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In general the transport equation for the Reynolds stresses (ujuj) can be written as; 

Qj = Dij +Pij +Fij - £ij + Rij ( 1 4 ) 

where C,y, Dy, Pij, Fij and £ represent convection, diffusion, production, pressure-strain and 
dissipation terms. The term Rij results from the transformation of the equation from plane to 
axially symmetric conditions and swirl. In Cartesian coordinates, the above terms are summarized 
below for each stress component; 


• u 2 - equation 


Cn = 7^ (prUu2)+ ~r^ (prVu2) 

_ Id , _ —5- k du 2 „ — k du 2 7 

z>;;= - r - 3 x lP '- C k ui — S ; + prC t uv— S rl 


1 d r ^ — k du 2 . _ ^ —k du 2 , 

+ -rS-lPrCtuv-^ + prC^-^l 

— dU—dU. 

Pyy= -2r(u 2 -^ + uv- w ) 


<Pll= -pC!j(u 2 -- 3 k)-C 2 (P n -- 3 P) 

+ P C lw \ [ - 2 U 2 fx + V 2 fy - uifxy ] 

+ C 2w [ 2 C 2 (P]] - | P)f x - C 2 (P 2 2 - | P)fy + C 2 P nfxy ] 

2 

£11= - jp£ 


• v 2 - equation; 

C 2 2= \.^(prUv 2 ) + - r ^(prVv 2 ) 
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P>22 = 


Id , „ k — dv* — dv' , , 1 d r . „ K ,— dv* . - 

73^ P rCk ~ (u ~d^ + uv ~dF )1 + rdr [ P rCk e (uV dx + V 

o <?V -y<?V — w . 

- 2 P( uv ~df + v If - VW T } 

-C lP l(v2-^k)-C 2 (P22-]P) 

+ p C] W ^( U 2 f X - 2 V 2 fy -uv fxy ) 

+ c 2w [ -C 2 (Pll - f Pj/x + 2 C 2 fP 2 2 - f P)fy+C 2 Pnfxy 

k . — ,o, ~ d , _ „ k — vw 

—1~ - T7E( P^k — 

e 


k y dv 2 


■ dv 2 


1 d 


k , — <?v 2 


P 22 = 
^22 : 


#22 = 


MW j 

r ' 


2 C k p^-\i<P c k ±(™?)- 2 i<pC k 

. 2 pCk ^. 2 pCk ^- 2 P C t ^ + 2 P ™ 


equation 

C 33 = L r ji(PrU^)* L r ^( P rV^) 


D 33 = ~ 


Id r k —dw 2 — dw 2 

-r7 i lprCk e <u ~^ 


“ v ~ 3 r, l 

dw 2 , 


1 d r k f dw 2 7 vrr x f 

+ 7 jp*>'-CW«v-j E - + vi sr * 1 


dW — dW — 7 V . 

P 33 = -2 p(uw-^ + vw-^r + w 2 -) 

F 33 = - p C] | (w 2 - | k) -C 2 (P 33 - 2 jP) 

+ p C/w | / m 2 f x + v 2 f y + 2 UVfxy ] 

- C 2 C 2w [ (P n -\p)fx + (P22- f P)fy-2C 2 P 12 fxy ] 


R 33 = 2pc k -^w 2 ^ + Ijif P c k~ (vw) 2 ) + 2-^( pCk-^uw^y) 



•mv - 


• vw 


, _ k (w 2 ) 2 - - kuwdvw „ „ fci — dvw _ — W 

- 2 p Q -- L 72 + 2 P Ck + 2 P Ck £ r VW dr ' 2p VW r 


equation 


Cl 2 = 
D]2 = 


- r ^(prU uv) + ~ r J; (prV mvJ 


7 <9 


it , — r <9wv — duv 


1 d 


uv 


l U r ^ o UL * v v \ i . * v r ~ r' / 

-rTx‘P rCt - (U ~ST +UV ~Si r,1+ ~r3r lprCk ~ (U ''Tx 


+ V z 


duv 

dr 


)] 


„ ,->dV —V -^dJJ —w, 

P 12 = - P ( u : 2 fa - UV - + v ; 2 - UW — ) 

£ 

012 = - pCij^uv - C 2 P 12 

- | pC Iw | [ UV(fx +fy) + (U 2 + V 2 )fxy ] 

+ | C 2 C 2w [ (P]l + P22 - 1 P)fxy + Pl2 (fx +fy) ] 


_ _ k uv 

Rj2 = - P Ck - -J - 
£ r 


d . _ k ( uw ) 2 . 7 <9 „ k . 

Tx< P Ct Z~r > - -r^<P C l-™™> 


_ 7 £ , — duw 

-pc k --(uw- 


_ dw 4 —W 
+ vw -gjr )+ puw — 


- equation 


C 23 = 
D 2 3 = 


p r U vw) + ^( p r V vw ) 

ld < -.*/ 3 rivw + — 


L rf x [p rCk ~ (u 2 ^r + 


dr 


1 d 


k , — <9 


+ 7 drtP r %( uv ~^ 


VW 


4* 


VW 


_<?W . — 

dr 


p 23 = ■Pl uv ^ + uw 7 2 


+ V z 


+ VW 


<9V 


)] 


V —,w 


dr 


+ vw — - w^ 
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£ 3 £ 

023= - pC!^vw-C2P23- 2 P C ^ k < uw fxy + vw fy) 

+ | C 2 C 2w (Pl3fxy + P23fy. ) 

— — , W ^ f — d — d ^ k — (v 2 -w 2 ; 4 

/? 2 i = -P (v 2 - w 2 ) — + p Cic ~ ~ vw (v 2 - w 2 ) + fa(p Cic~ uw ) 

£ £ 

+ jjfr(P C k ~vw(v 2 -w 2 ))-4 pC k -~ r vw^y + pCk^-uw-^fv 2 -w 2 ) 
£ £ £ 


♦ww-equation 


Cl3 - 
D]3 = 


~rlx ( ' pr U UW)+ ijfc(P rVuw ) 

1 d , „ k — dim . — duw 4 7 

- r Sx lprCk i <u ~sr* m ~s r)1 


1 d r _ k , — duw — 7 duw . , 

+ -r* [ P rCk ~ (uv -zr + v ~dT )] 


P 13 = 
P 13 = 


dW —dW 


, . v ,r dU UV . 

-P^^ + w "ZF + vw ‘SF' l<w '5 r; 

£ 3 £ 

-p Cj^uw- C 2 P 13 ~ 2 P Clw l uwf x 

3 £ 3 3 

- 2 P Clw J- vwfxy + 2 c 2w Q Pl3fx + 2 C 2 C 2 W P23 fxy 

— W „ k 1 — duv d , „ k — uv 

-puv T + pC k --vw^r+ Tx ( pCk-uw- ) 

Id , „ k , „ kuw duv _ „ k — uw 

+ 7 dr { P Ck Z VWUv) + pCk ~~r^- pCk ~ W 


dV 


Rl3 = 


The turbulence energy dissipation rate e is determined from its own transport equation; 


1 dprUkE 1 d 
r dx k 


ds 

r^/ rC &l'W‘'fo, , + C *i Pk - C * P T 


( 14 ) 


where the constants C £} and C £] have values of 1.44 and 1.92 respectively. 
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The terms f x , f y and f xy appearing in the stress-equation are tied to the orientation of the wall 
through the wall-damping function / and will be explained later in the wall reflection treatment 
section. 


5.3 Boundary Conditions 


To solve the transport equations for the Reynolds stresses, boundary conditions for the stresses are 
needed. In the present module the log-law based relations are used to bridge the gap between the 
fully turbulent and viscous near-wall regions. Boundary values for the stresses can be derived by 
applying the Reynolds stress equations to the near-wall equilibrium flow. It can be shown that the 


stresses are related to the turbulent kinetic energy = C,-,- k, where C,. are constants to be 
determined. Consider as an example, the log-layer turbulent flow, where S = , where u r 


is the friction velocity and K is Von Karman constant. In the log-layer, the limiting form of the 

stress equation is obtained by neglecting the convective terms and equating the production to 

dissipation and setting the wall-distance function /= 1, hence the molecular and turbulent diffusion 

terms can be neglected. Consequently, the normal stress equation for the wall-normal component 

2 

when simplified with &22 = Pj £ 1S ; 


^ 2 (-1 + Ci+ C 2 - 2C 2 C 2w ) r 

k ~ 3 (Ci + 2C lw ) 022 


(15) 


From experimental data, Lien & Leschziner [4] reported a value of C 22 ~ 0.249 for near wall 
equilibrium turbulence. The most frequently used value of Cj = 1.8 and C 2 = 0.6, and from 
Gibson and Launder [2] Cj w = 0.5 and C 2w = 0.3. Substituting these values into equation (13) 
give a value C 22 = 0.247 which is close to the experimental value. Similarly, these constants also 
give C]j = 1.09, C33 = 0.654 and C/ 2 = -0.255. 
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5.4 Numerical Procedure 


The conservation equations for the Reynolds stresses and the energy dissipation are integrated over 
control volumes after transformation of the Cartesian form to body-fitted no-orthogonal 
coordinates. The equation governing the transport of a scalar property 0, which stands for the 
Reynolds stress components and the energy dissipation equation can be written as; 

(pu ™ 0 - w ft J = s* ( 16 ) 

where represents the curvilinear coordinate frame and J is the Jacobian of the coordinate 
transformation, and represents its cofactors and qm represents the diffusion flux. Equation (16) 

is then integrated over discrete control volumes where the dependent variables on the volume faces 
are approximated by finite-difference representation. 

In general the diffusion term is represented as 

a?) 

where Tip is the diffusion coefficient. 


The tensorial form of the diffusivity due to Daly and Harlow [3] is adopted as; 


r<p — p r Cj Um u l 

£ 


(18) 


instead of the isotropic diffusivity (T# = fl t /(J<p )• Utilizing the equilibrium assumption and 
experimental near-wall stress data, the constant C s is taken to be 0.22 for the Reynolds stress 
equations and 0.18 for the turbulent energy dissipation equation. The diffusion term is discretized 
with a second-order central differencing scheme, while the convective terms are discretized using 
first or second order upwind differencing scheme. 


A special discretization practice for the Reynolds stress gradients is introduced into the finite 
volume procedure with colocated storage arrangement. This is necessary to avoid the problem of 
mean velocity-Reynolds stress decoupling that can lead to oscillatory solutions or even divergence 
of the iterative solution algorithm. The procedure adopted in the present work differs from that of 
Obi & Peric [12] and that of Lien and Leschziner [4] by accounting for all the driving forces of the 
Reynolds stresses and not only those given by the gradient-diffusion type process. To illustrate the 
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origin of the problem, consider the Reynolds stress gradient terms in the axial momentum equation 
in 2D Cartesian uniform grid for simplicity 

du^ Juv 

' dx ' By 

Now integrating over a control volume surrounding node P (cf. Figure 1 ) yields; 

" J l~^x~ dxd y-\ dxdy = -[(u 2 e -u 2 w ) Ay p + ( uv n - uv s ) Ax p ] 

Now if the cell face values of the u 2 e , u 2 w , uv n and uv s are evaluated with linear interpolation, the 
stress difference expression become 

_ l( 'i £ j^w )Ayp + ( uv N ^s )Axpl 

and since no P -node shear stress appear in the resulting expression, a chequer-board oscillation, 
similar to that played by the pressure field appear. Therefore, a non-linear interpolation scheme is 
needed to avoid these odd-even oscillations in the same context of Rhie and Chow [13] for cell face 
velocities. This means that any cell-face velocity is not merely sensitized to the pressure differences 
centered on that face but also to the Reynolds stress differences. 

Consider the descretized equation for the axial normal stress component u 2 in general non- 
orthogonal coordinates; 

A p u 2 p — E_^Ai u 2 /+ SiB (19) 

where n stands for the cells E,W,N and S neighboring P, A, are the coefficients for the 
neighboring cells and Su 2 is the source term that includes production, dissipation and pressure- 
strain redistribution terms as; 

S? -Pll - § pe+ 

where j combines Rotta's stress isotropization model and isotropization of production model 
and related wall-correction terms due to Gibson and Launder [2]. 

011= -pCij c (u 2 -^k)-C 2 (Pn-^P) 

+P C ]w j ( - 2u 2 f x + V 2 fy - UVf X y ) 

+ 2 C 2 C 2w ( Pll -jjP )fx - c 2 C 2w (P 22 - f P)fy + C 2 C 2w Pnfxy (20) 


Rearranging the production terms that contribute to the stress generation and noting that 
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( 21 ) 


P= ^Pkk, then; 

Su 2 = APu + BP 22 + CP 33 + DP 12 + $77 

where 

A = / - J C 2 + ^ C 2 C 2 wfy + J C2 C 2 wfx 
B = 1 - 3 C 2 + j C 2 C2wfy "t 3 C 2 C2wfx 
C = jC 2 + - 3 C 2 C 2w (fy ' 2f x ) 

D = C 2 C 2 wfxy 

and S]j contains the remaining terms. 

Substituting for the production terms Pjj, P22 > P33 ar >d ^ 72 . then equation ( 19 ) becomes; 


u 2 P = Hp + 2pA [ u 2 (D ] AUt + D 2 AU r l) + uv(EiAU'l + E2AUS)]p 

_____ — e y yyV W 

+ 2pB [ uv ( Dj AV % + D2 AVI) + v 2 ( Ej AV % + E2 AV $ ) ]p 

_ , , k w 2 W , 

+ 2pC [ uw( Dj AWt + D2AW 7 ! ) + vw ( E] AWS + E 2 AW$ ) - — 


+ pD [ u 2 ( D 1 AV % + D2 AV r l) + v 2 ( E] AU 7 + E2 AU £ ) + r ] p 

. ^77 


( 22 ) 


where 

H P = L-A^i/Ap 

D 1 = - Ay^/Ap , D2 = /Ap, 

E] = - Axp/Ap and £ 2 = Axj,/ Ap 

here Ay^ = ( V« -y s ) > &y\ - (ye - yw ) » etc 
and At/ £ = ('t/p - Up), AV £ — (Vp - Vp ) , etc 
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Now, performing the interpolation practice to obtain east cell-face value of the normal stress ( u 2 e ) 
we obtain; 

u 2 e = <u 2 p> - < 2pAu 2 D i AU % > - < 2pAuv £2 &U ^ > 

- < 2pBuv Dj AV Z > - < 2pBv 2 E2 AV £ > 

- < 2pCuwD] AW > - < 2pCvwE2 AW £ > 

- < pDu 2 D] AV > - < pD v 2 E 2 AU % > 

+ < 2pAu 2 D] > AU £ + < 2pA uvE 2 > AU % 

+ 2pB wv D 1 > AV % + < 2pV v 2 E 2 > AV % 

+ < 2pC uwDj > AW £ + < 2pCvw E 2 > AW % 

+ <2pD~u 2 Dj > AV$+ < pD~v2 E 2 > AU $ (23) 

The brackets < and > denote linear interpolation. For instance, on the east face 

Axp 

<<£>>= (1-M <Pp +/f &e where, ft 

Axp+Axp 

Similar expressions can be obtained for u 2 w „ u 2 n , u 2 s , uv w , uv n> and uv v which are then used to 
calculate the Reynolds stress gradients in the discretized axial momentum equation. Similarly, 

expressions for uv e , uv w , , uv n , uv s> ,v 2 e , v 2 w , v 2 s and v 2 n can be obtained for the stress 

gradients in the radial momentum equation and uw e , uw w , uw n , uw s , v w e , vw w , v w n and vw s 
expressions to evaluate stress gradients in the azimuthal momentum equation. 

5.4.1 Wall Reflection Treatment 

The wall reflection terms <P ij w appear in the pressure-strain term correlation as wall correction 
terms (0^- ]w and 0y 2w ) to counteract the tendency of 0^/ and <P ij2 to isotropise the 
turbulence structure. Special consideration is given to the wall proximity effects on the 
redistribution process <P ( j w with relation to the local orthogonal coordinate system at the wall, cf. 

figure 2. 

At a wall, turbulence approach a state of strong anisotropy associated with the tendency towards a 
2D turbulence. The wall-reflection terms ensure that normal stress normal to the wall is not too 
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high. For body-fitted coordinates, there is a need to consider the tensorial form of the wall 
reflection terms since they are tied to the orientation of the wall through the damping function term 
(eq. 12). For a curved surface (figure 2), the wall normal vector n = rtjij + n ^2 > where ij and 1*2 
are unit vectors in Cartesian coordinates. The Cartesian components of the wall-distance function f 
are given as; 

fx = n )f fy = n2 2 f ™&fxy = n l n 2f 

where f x = n 2 (cP' 75 k 1 - 5 / K £) / L n for example and L n is the normal distance from the wall. 

The Reynolds stresses close to the wall are transformed from wall coordinates to Cartesian 
coordinates by appropriate vector decompositioning to give; 

u 2 = u 2 tj + v 2 rij + 2 uv t]Ti i 
v 2 = u 2 t 2 + v 2 n 2 + 2 uv t 2 ti 2 
w 2 — w 2 

uv = u 2 tj t 2 + v 2 njn 2 + uv ( tj n 2 + t 2 nj ) 
v w = uw t 2 + vw n 2 
uw = uw tj + vw nj 


where u 2 , v 2 .... are the Reynolds stresses in Cartesian coordinates and u 2 , v 2 are the Reynolds 

stresses in wall-coordinate, n j, n 2 are the Cartesian components of the normal vector component 
and tj, t 2 are the Cartesian components of the tangential vector component. 


5.5 Module Evaluation 

The RSM module was tested at Rocketdyne after interfacing with the CFD solver REACT and at 
the University of Alabama at Huntsville (UAH) using own solver (MAST). The first test was on 
fully developed channel flow with length to height ratio of 50 and a Reynolds number of 2xl0 5 
based on the channel height. A non-uniform mesh of 101x41 was used with clustering at the walls. 
Figure 3 shows the fully developed mean velocity profile across the channel. Figure 4 shows the 
normal Reynolds stress profiles across the channel and figure 5 shows the shear stress profile. 
Similar results were obtained when the module was interfaced and tested independently at UAH 
using the MAST code. 
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The next test problem is that of the backward facing step of Driver and Seegmiller [14], The 
calculations were performed using a 101x41 grid points with clustering near the walls. The 
computational domain had a length of 50H (H is the step height) and a width of 9H. The 
experimental data were used to specify the inflow conditions for a channel flow calculation where 
the fully developed profiles at the channel exit were used as the inlet conditions for the backward 
facing step calculations. Fully developed flow conditions were imposed at the outflow boundary. 
The boundary conditions for the Reynolds stress equations were arrived at by using the log-law of 
the wall and assuming local equilibrium conditions close to the wall. It can be shown that the 
Reynolds stresses are related to the turbulent kinetic energy by; 

vijUj = Cy k (24) 

were Cy is constant. The Reynolds stresses at the vicinity of the wall used are 

IP = 1.098 k, v 7 = 0.247 k, uv = - 0.255 k and w 2 = 2k -IP - P — 0.654 k. 

Figure 6 shows the stream lines using Launder, Reece and Rodi’s model. The computed 
reattachment length was about 5.8H which is closer to the experimental value of 6.1H than the 
standard k-e model (5.35H). The figure also shows a small (turbulence driven) secondary flow 
region at the base corner of the step which cannot be predicted using the isotropic eddy-viscosity 
k-e model. Also, a smaller recirculation region is noted at the top lip of the step which is also 
driven by turbulence anisotropy (more refined grid may be needed to isolate and study this region). 
Figure 7, shows the mean velocity profile across the channel at four step heights downstream of 
the step as compared with the standard k-e turbulence model predictions. The axial normal 

turbulent intensity (IP ) profile across the channel at x/H=4 is shown on figure 8 and the radial 

normal turbulent intensity (P ) is shown on figure 9. The shear stress (uv ) profile across the 
channel at x/H=4 is also shown on figure 10. The results show that the module predicts improved 
results using the RSM model as compared with the standard k-e model. 
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Figure 3. Mean velocity profile across the channel 



Figure 4. Turbulent intensity profiles across the channel 
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Figure 5. Shear stress profile across the channel 











Figure 8. Turbulent intensity uu profile 

(normalized with U 2 ref ) at X/H = 4 
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Figure 10. Turbulent shear stress uv profile 
(normalized with ) at XIH = 4 



APPENDIX D 


2D/Axisymmetric Reynolds Stress Module Deck 

The 2D/axisymmetric Reynolds stress module is a FORTRAN source code to solve 
2D/Axisymmetric turbulent flow using the full Reynolds stress model based on Launder, Reece 
and Rodi[l] when interfaced with a main flow solver. The module consists of the main routine 
RSMOD that calls a number of subroutines to perform different functions that will be explained 
below. 

3.1 Subroutine RSMOD 

This is basically the main routine that reads through its argument list different variables from the 
calling flow solver which are described below. 

List of Argument Variable Names 

X Grid node locations in the x or ^-direction, dimensioned to X(NX*NY) 

Y Grid node locations in the y or T|-direction, dimensioned to Y(NX*NY) 

FX Interpolation factor in the x or ^-direction. 

FY Interpolation factor in the y or T] -directi on. 

ARE Control cell areas 

VOL Control cell volumes. 

R Radial distance in the axisymmetric geometry or 1 . for planar geometry. 

DNS Normal distance of a cell from the south-boundary dimensioned to NX. 

DNN Normal distance of a cell from the north-boundary dimensioned to NX. 

DNE Normal distance of a cell from the east-boundary dimensioned to NY. 

DNW Normal distance of a cell from the west-boundary dimensioned to NY. 

U Axial or ^-direction velocity, dimensioned to NX*NY. 

V Radial or vjr-direction velocity, dimensioned to NX*NY. 

W Tangential or azimuthal velocity, dimensioned to NX*NY. 

TE Turbulent kinetic energy, dimensioned to NX*NY. 

ED Turbulent energy dissipation rate, dimensioned to NX*NY. 

DEN Density (assumed constant for incompressible flows). 

FI Mass flux at cell faces in the x or ^-direction, dimensioned to NX*NY. 

F2 Mass flux at cell faces in the y or redirection, dimensioned to NX*NY. 
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VISCOUS 

VIS 

RESOR 

ITBS 


ITBN 

JTBE 

JTBW 

ITER 

ICAL 

AKSI 

RESTART 


Laminar viscosity. 

Eddy viscosity, dimensioned to NX*NY. 

Residual error for the equations solver, dimensioned to 8. 

Boundary condition flag along the south boundary dimensioned to NX and must 
have one for each boundary node set to: 1 -inlet, 2-outlet, 3-symmetiy and 4-wall 
e.g., for a wall boundary condition along the south boundary set ITBS to NX *4. 
Similarly for the other boundaries. 

Boundary condition flag along the north boundary, dimensioned to NX. 
Boundary condition flag along the east-boundary dimensioned to NY. 

Boundary condition flag along the west-boundary dimensioned to NY. 

Iteration number. 

= 1 for swirl velocity calculations, 0 otherwise. 

= l for axisymmetric flow, 0 otherwise. 

= l if calculations are restarted from a previous run, 0 otherwise. 


RSMOD starts by reading the turbulent flow constants, under-relaxation factors and 
Prandtl/Schmidt numbers for the k and £ equations. These are; 


CD1, CD2 

CMU, ELOG 
and CAPPA 
GKE 

ALFAKJE 

URFVIS 

SORKE(l-8) 

URFKEU-8) 

PRTKEO-8) 


Cj, C 2 


constants in the k and e -equations and are usually set to 1.44 and 1.92 
respectively. 

also constants in the k and £ -equations and are usually set to 0.09, 9.8 and 
0.42 respectively. 

is set to 1 for second-order upwinding of the convective terms in the 
transport equations. 

is the iteration parameter used in the k and £ -equation solver, 
is the underrelaxation factor of the viscosity near the wall. 

are the degree of accuracy for the k, e, u 2 , v 2 , w 2 , uv, vw, and uw- 
equations solver respectively. 

are the underrelaxation factors for the k, £, u 2 , v 2 , w 2 , uv, vw, and ww- 
equations respectively. 

are ratio of Prandtl to Schmidt numbers used in the k, £, u 2 , 

v 2 , w 2 , uv, vw, and uw -equations respectively, 
are constants in the RSM model. 
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Cip and C 2 P 


are the two constants in the wall-reflection terms of the pressure-strain 
redistribution term. 

Ck and C e constants in the diffusion term of the k and s -equations. 

CUU, CVV, CWW, are the constants multiplying the kinetic energy for the stress values near 

CUV, CVW, CUW the wall. 

WREFON = 1 if the wall reflection terms of the pressure-strain term are to be included, 

= 0 otherwise. 

All variable dimensions considered are one-dimensional. The position of any node is defined as IJ 
= (I,J) = (I-1)*NJ +J, where NI and NJ are the number of grid nodes in the X and Y-directions 
respectively. It is assumed that grid related data such as cell areas, volumes and interpolation 
factors be passed to the module from an external grid generator. 


Subroutine CALPIJ 

This subroutine calculates the production terms of the individual stress components. 

Subroutine CALUIUJ 

This subroutine solves the transport equations for the turbulent energy (IPHI=1), energy 
dissipation. (IPHI=2) and Reunolds stresses (IPHI=3, 4, 5, 6, 7, 8 for 

u. 2 ,v 2 ,w 2 , ~uv, vw, and uw ). Daly and Harlow [3] gradient stress diffusion form is used in the 
module instead of the simplified isotropic diffusivity form. This subroutine calls MODUIUJ 
subroutine that sets the appropriate boundary conditions for the Reynolds stresses. The set of 
algebraic difference equations are then solved using Stone's strongly implicit solver SOLSIP. 

Subroutine MODPIJ 

This subroutine modifies the production terms near the wall using the near wall region model. 

Subroutine MODUIUJ 

This subroutine calculates the near wall boundary conditions for all the variables. 
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Subroutine SOLSIP 


This subroutine solves the system of linear algebraic equations for all the variables using Stone's 
Implicit Procedure. 

Subroutine WALREF 

This subroutine calculates the wall reflection terms in the pressure-strain redistribution correlation. 
It calculates the wall unit normal vectors and the normal distance away from the wall. This is 
needed to resolve the wall tangential and normal velocity components that are needed to obtain the 
near-wall values of the Reynolds stresses. 

SUBROUTINE WALPARA 

This subroutine calculates the normal and tangential wall unit vectors. 
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6.1 Introduction 


In this section a description of the standard k-e turbulence model that is coded as a self contained 
computer program to compute turbulent flow quantities in three-dimensional, body-fitted geometry 
is given. Module structure and variables used are given in the Appendix. The module was 
successfully tested as a self-contained unit using the REACT code[l]. 


6.2 Theory and Model Equations 


The k-e turbulence model used is based on the standard two equation k-e model of Launder and 
Splading [2], For a steady, incompressible flow the transport equations for the turbulent kinetic 
energy k and energy dissipation £can be written in generalized Cartesian coordinates as; 



L A 

r dxj 


(H+ )+G-pe 

°k J 


( 1 ) 


he Id. Pt de 
t=~r + 


( 2 ) 


where G denotes the rate of production of turbulent kinetic energy and is expressed as; 


„ dU, , dU; 



) 


The empirical constants, a <J E , Cj and Cj have values 1.0, 1.0, 1.44 and 1.92 respectively. 


The above equations are valid only in the fully turbulent region away from the wall. Therefore the 
wall function method (similar to that described in Chapter 2 for the 2D k-e module) is used to 
model the damping effects of the thin sublayer region close to the wall. 



6.3 Module Evaluation 


The 3D k-e turbulent module was evaluated by interfacing the module with the REACT code as the 
CFD solver and producing the same results that were generated previously with the full REACT 
code for a centrifuged impeller calculations (Chen et. al [3]). Figure 1 shows the grid topology of 
the impeller studied with the shroud removed and Figure 2, shows the reduced pressure plot . In 
general Chen et al's calculations showed good comparisons with experimental data obtained from 
laser velocimetry in a water test rig. 
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Figure 1. Impeller grid topology 




Figure 2. Reduced pressures 


APPENDIX E 


3D k-e Turbulence Module Deck 

This module consists of two separate programs KEMOD3 and MODIFY, which have to be linked 
to the main flow solver. A description of each file will be given next. 


Program KEMOD3 

This is basically the solver for the k and e - transport equations. It reads through its argument list 
different variables from the calling flow solver. These variables are described below. 


List of Argument Variable Names 

NIM Number of cell nodes in the I- or ^-coordinate lines, (input from the flow 

solver) 

NJM Number of cell nodes in the J- or ri-coordinate lines, (input from the flow 

solver) 

NJM Number of cell nodes in the k- or ^-coordinate lines, (input from the flow 

solver) 

X Grid node locations in the x or ^-direction, dimensioned to X(JXYZ) 

(JXYZ=NX*NY*NZ) (input from flow solver) 

Y Grid node locations in the y or r| -direction, dimensioned to Y(JXYZ) (input 
from flow solver) 

Y Grid node locations in the z or ^-direction, dimensioned to Y(JXYZ) (input 
from flow solver) 

U x-direction velocity (u), dimensioned as U(JXYZ) (input from flow solver) 

Y y-direction velocity (v), also dimensional as V(JXYZ) (input from flow solver) 

W z-direction velocity (w), dimensional W(JXYZ) (input from flow solver) 

TE Turbulence kinetic energy k , dimensioned TE(JXYZ) (calculated in the module 

and returned to the flow solver) 
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ED 

URFK 

URFE 

PRTK 

PRTE 

G 


FI 

F2 

F3 

ITER 

VISCOS 

VIS 

Cl 

C2 

CMU 

BCFE 


BCFW 

BCFS 


Turbulent energy dissipation rate e, dimensioned ED(JXYZ) (calculated in the 
module and returned to the flow solver) 

Under-relaxation factor for k -equation (input from flow solver) 
Under-relaxation factor for e-equation (input from flow solver) 

Prandtl/Schmidt number for turbulent energy-equation, assumed known (input 
from flow solver) 

Prandtl/Schmidt number for turbulent energy dissipation equation, assumed 
known (input from flow solver) 

= 1 .0 if second order upwinding is desired 
= 0.0 if first order upwinding is used 

(input from flow solver. Usually calculation of k and e are not very sensitive to 
the order of upwinding used) 

Mass flux variable at cell faces in x- or ^-direction, dimensioned Fl(JXYZ) 
(input from flow solver) 

Mass flux variable at cell faces in y or q -direction, dimensioned F2(JXYZ) 
(input from flow solver) 

Mass flux variable at cell faces in z or ^-direction, dimensioned F3(JXYZ) 
(input from flow solver) 

Iteration number (input from flow solver) 

Dynamic viscosity (input from flow solver) 

Eddy viscosity, dimensioned VIS(JXYZ) (calculated in the module and returned 
to the main solver) 

Turbulence model constant, Ci (input from flow solver) 

Turbulence model constant, C 2 (input from flow solver) 

Turbulence model constant, (input from flow solver) 

Boundary condition flag along east boundary (or y-z plane). It must have one 
for each boundary node set to: 1-inlet, 2-outlet, 3-symmetry and 4-wall e.g., 
for an outlet boundary condition on the east boundary set EBCE to (NY*NZ)*2, 
and similarly for other boundaries, dimensioned BCFE(JYZ=NY*NZ) (input 
from flow solver) 

Boundary condition flag along west boundary, dimensioned BCFW(JYZ) 
(input from flow solver) 

Boundary condition flag along the south boundary, dimensioned 
BCFS(JXZ=NX*NZ) (input from flow solver) 
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BCFN Boundary condition flag along north boundary, dimensioned BCFN(JXZ) 

(input from flow solver) 

BCFB Boundary condition flag along bottom boundary (or x-y plane), dimensioned 

BCFB(JXY=NX*NY) (input from flow solver) 

BCFT Boundary condition flag along top boundary (or x-y plane), dimensioned 

BCFT(JXY=NX*NY) (input from flow solver) 


The module is interfaced with the main flow solver by a call to KEMOD3 with its arguments. For 
iterative flow solvers the module is called within the iteration sequence after the solution of the 
momentum equations where the mean velocities are passed to the module. There are different flow 
solvers utilizing different schemes from staggered to nonstaggered grid arrangement and for 
nonorthogonal coordinate system there are at least three alternatives to the choice of the velocity 
components; 

i. Cartesian velocity components 

ii. Contravariant velocity components 

iii. Covariant velocity components 

The Cartesian velocity components are the most widely used and have the advantage of simple 
formulation of the governing equations. Whatever the arrangement used, mass fluxes at cell faces 
are required and passed to the module as FI, F2 and F3 in all directions. The location of other 
variables such as k and e are at the cell center or cell nodes. 

The module starts by reassigning variable names passed to it from flow solver to names that are 
shared with the different subroutines of the module in a common statement file "kemod.h". Then a 
check is made if it is the first iteration in which case the grid file "GRIDG" is called -after passing 
the grid node locations X, Y and Z- in order to calculate grid related quantities which will be 
explained later. The need to call GRIDG can be waived if all the grid data are passed to the 
module. That is all the information about the grid such as interpolation factors FX, FY and FZ, cell 
volumes (VOL) and normal distances of first grid point from grid boundaries (DNS from south 
boundary, DNN - from north boundary, DNW - from west boundary, DNE - from east boundary, 
DNB - from bottom boundary and DNT - from top boundary). 

After this a call to subroutine CALCE is made to calculate the turbulent kinetic energy k (with the 
identifier IPHI=1). The energy dissipation equation is solved next by a call to subroutine CALCE 
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again with the identifier IPHI=2. The turbulent viscosity is updated next by calling subroutine 
MODVIS. A brief description of each subroutine is given next. 


Subroutine GRIDG 

Before calling this subroutine, the coordinates of all grid nodes, defined in reference to a fixed 
Cartesian coordinate frame are read. Figure 3 shows the position of cell and grid nodes. The west- 
to-east, south-to-north and bottom-to-top directions correspond to the ascending indexing order of 
i, j and k, respectively, forming a right-handed coordinate system. 

This subroutine is called only once to calculate coordinates of grid nodes (intersection of grid lines) 
and geometrical properties of the grid (cell volumes, interpolation factors, normal distances of 
near-boundary cell nodes from boundary). All variables including grid node coordinates are 
converted to one-dimensional arrays. The position of any node in one-dimensional array is 
therefore defined as; 

IJK = (I-1)*NJ + (K-1)*NI*NJ + J 

where NI, NJ and NK are the maximum number of grid nodes in the i, j and k directions 
respectively. 

The actual number of grid nodes is one row and one column less than for all cell nodes. For I = 
NI, J = NJ and K = NK fictitious grid nodes are introduced which have the same coordinates as 
actual nodes on NI-1 in I-direction, NJ-1 in the J-direction and NK-1 in K-direction. 

The subroutine then calculates interpolation factors which are associated with cell nodes and are 
used in the main program to calculate values of dependent variables at locations other than cell 
nodes (cell centers). Cell volumes are calculated next followed by calculations of normal distances 
of near-boundary nodes from all the six outer boundaries. 
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Subroutine CALCE IPHI. IPHI) 


This subroutine solves the linearized and discretized transport equations for the turbulent energy k 
and the energy dissipation rate e. The two dummy parameters in the calling statement, PHI and 
IPHI, represent arrays containing dependent variables for which the equation is to be solved. The 
subroutine sets up the convective and diffusive coefficients over the entire field, then it calculates 
the source terms for either k or £ transport equations. A call is made to MODKE or MODED in 
order to modify the sources for k and £ equations respectively. 

The discretized equations have the form 


A p <P p = X A i 0 i + 

i=EWNSTB 

where the coefficients A, (i=E,W,N,S,T,B) contain both the convective and diffusive fluxes, these 
equations are assembled and solved by calling subroutine SOLSIP which is based on Stone's 
Strongly Implicit Solver [4], 

Subroutine SOLSIP 

This subroutine solves the system of linear algebraic equations for k and £ using Stone's Implicit 
Procedure [4]. The array RES (IJK) is used to store the residuals. The sum of absolute residuals 
"RESl" calculated in the first pass through this part of the routine is used as a measure of 
convergence of the solution process as a whole and this value is stored in RESOR (IPHI). This 
variable RESOR (IPHI) is passed to the main solver and if desired can be normalized and 
compared with the maximum error allowed there. If necessary inner iterations counter L and the 
sum of absolute residuals RES 1 are printed out to monitor the rate of convergence of k and £ 
solution. If the ratio RSM is greater than the maximum allowed for the variable in question, SOR 
(IPHI), and the number of inner iterations is smaller than a prescribed maximum, NSWP (IPHI), 
then the routine repeats the sequence of calculating the residuals, increment vectors and updating 
the dependent variable. 
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Subroutine MODVIS 


This section calculates the effective viscosity and is called after calculating k and e. At locations 
where e is close to zero (i.e., < 10' 30 ) viscosity is set to zero. A provision is made for under 
relaxing changes in effective viscosity which may help to stabilize oscillations and improve 
convergence rate. 

Subroutines MODK and MODED 

These subroutines are called from subroutine CALCE and they set the boundary conditions for 
k and £ . For the kinetic energy equation for example, the bottom boundary is checked 

first for one of the options below; 

(1) An inflow boundary BCFB(IJ) = 1 (IJ = (I-1)*NJ+J), where the source term is set to 
accept the inlet values at the x-y plane (bottom boundary K=l). 

(2) Outflow boundary BCFB(IJ) = 2, where zero gradient in the z-direction is employed. 

(3) Symmetry boundary, BCFB(IJ) = 3, where gradients normal to symmetry x-y plane are 
zero. 

(4) Wall boundary, BCFB(IJ) = 4, where the turbulent kinetic energy production (per unit 
volume) term GENTB(I) calculated form subroutine WALLFN in program MODIFY is 
added to the rest of the source term SU(IJK). 

Boundary conditions for the £ -equation are similar to those of k except at the wall where they are 
set to appropriate values for each near wall treatment. 

Program MODIFY 

This program is compiled separately and is called from the u, v and w momentum solver .It 
basically updates the flux source term of the discretized momentum equation due to wall shear 
stresses. If the u-momentum equation for example is discretized in the form 
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A p u p 


y, A( ui 

i=EWNSTB 


* 

+ 5 

u 


* 

where P, E , Wi N, S, T y B are cell nodes, and A and A/'s contain convective and diffusive 

coefficients. S u is the source term containing pressure gradients and cross-derivative diffusion 

terms and convective terms for second-order upwinding scheme. This source term is usually 
* 

linearized as S u = S u - B p u p The term B p is usually moved to the left hand side of the equation and 

* 

modifies the diagonal coefficient A p =A+B p9 and the equation can be written as 


Ap Un — ^ Ai Ui "t" Su 

y i=EWNSTB 

Then S u and B p are passed to subroutine MODIFY where they are modified if a wall is present 
(e.g., BCFB(IJ) = 4 for bottom boundary). 
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7.1 Introduction 


In this section a description is given of the three-dimensional Algebraic Stress turbulence Model 
(ASM) based on the work of Rodi [1]. The model is coded as a self contained computer program 
to compute turbulent flow quantities when interfaced with a CFD solver. Detailed description of the 
module structure, variables used and how to interface the module with CFD flow solvers are given 
in the Appendix. 

The module uses as input the mean flow properties, as computed by conventional CFD solvers, 
and calculates the Reynolds stresses, turbulent kinetic energy and the energy dissipation. It is 
structured to be self-contained and compatible with many CFD codes. The module has not been 
tested thoroughly due to the ending of the contract earlier than scheduled. Some testing of the 
module has been done at UAH but that also has been put on hold. However, the module as 
assembled is capable of interfacing with a number CFD solvers. 

The module computes turbulent flow quantities in three-dimensional body-fitted geometry with or 
without rotation about any one of the three axis. The standard wall functions is used for the near 
wall treatment. 

7.2 Theory and Model Equations 

The Algebraic Stress (ASM) module discussed here is based on the work of Rodi [1]. The idea is 
to simplify or truncate the Reynolds stress equation by approximating the convective and diffusive 

transport of the Reynolds stresses ujuj in terms of the corresponding transport of turbulent 
energy. This allows the transport equation for the stresses to be expressed as a set of algebraic 
formulae containing the turbulence energy and its rate of dissipation as unknowns in the form: 

k 2 

UiUj = JP-e) [ Plj ' 3 6ij £ + 0ij] 

where Pij = Production and P =T Pk k 

<Pjj = pressure-strain redistribution 

®ij = & ij, 1 + &ij, 2 + ®ij,1w + &ij,2w 
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Rotta's linear retum-to-isotropy concept for the non-linear part 

£ 2 

®ij,1 = -Ci ^ (ujuj - ~ k 6ij) 

is used and the "isotropization of production" concept for the linear "rapid" part 

®ij,2 = -C 2 ( Pij - ^ P Sjj) 

is used. Gibson and Launder [2] concept for the wall reflection terms is used as 
£ 3 3 

Oij.lw = Clw ( uku m nkn m 8ij - J riknj - j UkUj nkni )f 


*&ij,2w — C2w ( & km, 2 n k n m$ij ~ 2 ^0^,2 n k n j ~ 2 ®}k ,2 n k n i )f 


where (rip is the wall-normal unit vector in the i -direction. The wall-distance function (f) 

1 ( 3/2 

represents the ratio of the turbulence length scale (L e = ) and the wall distance and is given 


as 


< 75 k‘ s 1 

f= ( ) — 

J 1 K e ’ An 

with An being the wall-normal distance. 


The resulting set of algebraic equations for the Reynolds stresses can be arranged in the form 
Aij u 2 + Bij v 2 + Cij w 2 + Dij u v + E[j vw + Fy uw= Gij 

where A jj, Bjj, Cij, Djj, Ejj, Fjj, and Gij are functions of the mean and turbulent flow 
variables. 


The above equation can be solved iteratively in the main flow solver. However, the algebraic 
system of equations is stiff and convergence difficulties are encountered when solved iteratively. 
Therefore, the set of equations was cast in the general matrix form A T = B , where 
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3e 

d(J 

dv 

dv 


, dJ dv _ ^ 

dv dw 

du dw 

C 0 Oy 

+ 2 
2Xk 

dx 

dy 

dz 


l dy dx- 6CoQz 


2 ^-J7 + 6 

du 


3e dV 

dv 


dV dU ^ _ 

„ dV dW ^ 

dU 

dW 


dx 


+ 2 j 

2Xk dy 

dz 


2 dx~ dy +6CoQz 

‘ <?>- ' 6C °A 

- ( Tl 

+ dx ) 


du 


dV 

3e 

dw 

,dU dV t 

<9v ^ _ 

du 

dw 

6C 0 Q y 

dx 


-dy 

+ 2 
2Xk 

dz 

‘ dy + dx } 

** +2 ^7 +6C oA 

' dz* 

dx 

dV 0 


dU „„ 

0 


£ dU dV 

du ^ 

dV 



^ +2 

dy ' 2C °A 


A k + dx + dy 

dz + 2 c ° A 

dz 

- 2 C 0 Q X 

0 



dv 

Cq ^ 2 x 

dw ^ 

£ <?V oW 

dV 




dy + 2 C 0 A 

T z - 2 

d x - 2Co °y 

A* + 

dx* 

2 C 0 Q 


dw . 


0 

du 

’.C 0 £2y 

dW „ _ 

<*/ _ 

£ 

dU 

dW 

^' 2C «A 

Tz + < 

dy + 2C ° A 

^ + 2 QA 

Xk 

+ dx + 

!k 


T=[puu,pvv,pww,puv,pvw,puw ] 


B 


where X = 


pe 3 

X + 2(1 -C 2 ) 11,1w + ® 11,2w ) 

pe 3 

X + 2(1 -C 2 ) ( 022 ' 1w + °22,2w) 
pe 3 

X + 2(1 -C 2 ) ^ 33 ’ 1w + °33,2w) 

(1 -C 2 ) (® 12 ’ 1w + ®12,2w) 

JJTc^j @23,1 w + ®23,2w) 
(1-C 2 ) (® 13 '1w + &13,2w) 

I-C 2 


C-i-1 + 


pe 


The matrix was inverted at each iteration step to obtain a converged solution. 
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APPENDIX F 


3D Algebraic Stress Module Deck 

ASMOD is a FORTRAN source code to solve 2D/Axisymmetric turbulent flow quantities using the 
algebraic stress model when interfaced with a main flow solver. The module consists of the main 
routine ASMOD that calls a number of subroutines to perform different functions that will be 
explained below. 


Subroutine ASMMOD 


This is basically the main routine that reads through its argument list different variables from the 
calling flow solver which are described below. 

List of Argument Variable Names 


INTTASM 

NIM 

NJM 

NKM 

U 

LK 

FX 

FY 

FZ 

X 

Y 

Z 

VOL 

U 


Initialization parameter that writes and sets variables 
Number of grid nodes in the i (or x) direction 
Number of grid nodes in the j (or y) direction 
Number of grid nodes in the k (or z) direction 

LI(I)=(I-1)*NJ, dimensioned to NX. Calculated as in subroutine GRIDG of the 3D 
k-e module. 

LK(K)=(K-1)*NI*NJ dimensioned to NZ. Calculated as in subroutine GRIDG of 

the 3D k-e module. 

grid interpolation factor in the x-direction 

grid interpolation factor in the y-direction 

grid interpolation factor in the z-direction 

Grid node locations in the x or ^-direction, dimensioned to 

X(JXYZ=NX*NY*NZ) 

Grid node locations in the y or r)-direction, dimensioned to Y(JXYZ) 

Grid node locations in the z or ^-direction, dimensioned to Z(JXYZ) 

Control cell volume (similar to that calculated in GRIDG of k -e module) 
mean velocity in x or ^-direction, dimensioned to U(JXYZ) 

(input from the flow solver) 
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V 

w 

VIS 

TE 

ED 

U2 

V2 

W2 

uv 

vw 

uw 

GEN 

SUASM 

SVASM 

SWASM 

BCFW 


BCFE 

BCFS 

BCFN 

BCFB 

BCFT 

GENTW 

OMX 


Mean velocity in the y or T| -direction, dimensioned yo V(JXYZ) 

(input from the flow solver) 

Mean velocity in the z or ^-direction, dimensioned yo W(JXYZ) 

(input from the flow solver) 

Eddy viscosity 

Turbulent kinetic energy, dimensioned to TE(JXYZ) calculated in the module. 
Turbulent energy dissipation, dimensioned to ED(JXYZ) 

Normal Reynolds stress component u 2 , calculated in the module 
Normal Reynolds stress component v 2 , calculated in the module 
Normal Reynolds stress component w 2 , calculated in the module 
Shear stress component uv , calculated in the module 
Shear stress component vw calculated in the module 

Shear stress component uw , calculated in the module 
Turbulent energy generation term 

Source term for the U-momentum equation due to Reynolds stress gradients. 
Calculated in the module and passed to the main solver. 

Source term for the V-momentum equation due to Reynolds stress gradients. 
Calculated in the module and passed to the main solver. 

Source term for the W-momentum equation due to Reynolds stress gradients. 
Calculated in the module and passed to the main solver. 

Boundary condition flag along the west boundary (or y-z plane). It must have one 
for each boundary node set to; 1-inlet, 2-outlet, 3-symmetry and 4-wall. For 
example for an outlet flow condition on the west boundary set BCFW to 
(NY*NZ)*2, and similarly for the other boundaries, dimensioned to 
BCFW(JYZ=NY*NZ) (input from flow solver) 

Boundary condition flag for the east boundary dimensioned to BCFE(JYZ) 
Boundary condition flag for the south boundary dimensioned to BCFS(JXZ) 
Boundary condition flag for the north boundary dimensioned to BCFN(JYZ) 
Boundary condition flag for the bottom boundary dimensioned to BCFB(JXY) 
Boundary condition flag for the top boundary dimensioned to BCFT(JXY) 
Turbulent generation terms calculated from the wall functions close to the wall in 
the west direction. Similarly for the other GENTE, GENTS.... 

Frame rotation term in the x-direction. Similarly OMY & OMZ in the 
y and z-directions respectively 


-96- 



DENSIT Constant density 

VISCOS Kinematic viscosity 

All dimensions considered are one-dimensional. The position of any node is defined as 
IJK = (I,J,K) = (I-1)*NJ + (K-1)*NJ + J, where NI, NJ and NK are the number of grid nodes in 
the X, Y and Z-directions respectively. It is assumed that grid related data such as control volumes 
and interpolation factors be passed to the module from an external grid generator, similar to the one 
listed in the 3D k-e module (Chapter 6). 

Subroutine CALPIJ 

This subroutine calculates the production terms of the individual stress components. 

Subroutine CALUIUJ 

This subroutine calculates the individual stress component from its algebraic equation. It sets the 
coefficients of the algebraic stress equations which are solved implicitly at each iteration step by 
inverting a 6x6 matrix. 

Subroutine SORUVW 

This subroutine calculates the source terms needed in the momentum equation of the main CFD 
solver due to Reynolds stress gradients. 

Subroutine SOLV 

This subroutine is a Gaussian elimination solver to invert a 6x6 matrices. 

Subroutine WALSTRS 

This subroutine calculates the Reynolds stresses near the walls based on wall functions. 
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CHAPTER 8 


Related Publications and Presentations 

During the course of this work, some related turbulence modeling work was published or 
presented at different meetings. Copies of these papers are listed in this chapter, they include; 

(1) A. H. Hadid ad M. M. Sindir "Comparative study of advanced turbulence models for 
turbomachinery” NACA CP-3174, 1992. 

This paper was presented at the advanced Earth-to-Orbit propulsion technology conference 
held at NASA Marshall Space Flight Center in Huntsville, Alabama on May 19-21, 1992. The 
work tests different correction to the standard k-e turbulence models that accounts for streamline 
curvature and rotations using different near wall treatments. 


(2) A. H. Hadid, M. E. DeCroix and M. M. Sindir "Advanced turbulence models for 
turbomachinery" NASA CP-3221, 1993. 

This paper was presented at the eleventh workshop for compuatational fluid dynamics 
applications in rocket propulsion held at NASA Marshal Space Flight Center in Huntsville, 
Alabama on April 20-22, 1993. The paper outlined the progress of the 2D/axisymmetric single and 
multi-scale k-£ turbulence module deck developments. 


(3) A. Hadid, M. Sindir, C. Chen and H. Wei "Computations of confined swirling flows with 
high order turbulence models in a modular form" 

This paper was presented at the twelfth workshop for computational fluid dynamics 
applications in rocket propulsion held at NASA Marshal Space Flight Center in Huntsville, 
Alabama April-May 1994. The paper presented the status of the 2D/axisymmetric second order 
closure models using the algebraic and the full Reynolds stress models. 

(4) A. H. Hadid, M. M. Sindir and R. I. Issa "A numerical study of two-dimensional vortex 
shedding from rectangular cylinders" published in the CFD Journal July, 1992. 
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This paper presents a test for an anisotropic k-e turbulence model. This model is an 
improvement on the standard k-e model since it can predicts Reynolds stress anisotropies without 
the need to solve additional equations for the stresses. 


(5) A. H. Hadid, N. N. Mansour and O. Zeman "Single point modeling of rotating turbulent 
flows" Proceedings of the 1994 summer program, CTR, NASA Ames/Stanford 
University. 

This paper tests a new one-point closure model that incorporates the effects of rotation on the 
power-law decay exponent of the turbulent kinetic energy. A modification to the e -equation 
proposed by Zeman using large eddy simulation results was used. A new definition of the mean 
rotation was proposed based on critical point theory to generalize the effects of rotation on 
turbulence to arbitrary mean deformations. 
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COMPARATIVE STUDY OF ADVANCED TURBULENCE MODELS FOR TURBOMACHINERY 


A. H. Hadid and M. M. Sindir 
CFD Technology Center 
Rocketdyne Division, Rockwell International 
Mail Code IB39, 6633 Canoga Avenue 
Canoga Park, CA 91303 


ABSTRACT 

Development and assessment of the standard k— E turbulence model for rotating flows with different near- wall 
treatments is presented. These include the standard wall function^, Patel’s two-layer model (2 \ and Lam and Bremhorst< 3 > 
low-Reynolds number model. Two test cases were chosen to validate these models for rotating flows. The first, from Daily 
and Nece^ f is for a rotating disk cavity in which recirculation and secondary flows are induced by the rotating element The 
second case is that of a confined double concentric jets with a sudden expansion by Roback and Johnson^. 

It is shown that near-wall effects are important close to rotating walls and that the two-layer model behaves better 
than the other two near-wall models. For confined swirling flows with fixed walls, the near wall effects are of secondary 
importance to the Reynold’s stress anisotropy. 


INTRODUCTION 

Accurate predictions of turbulent flows are crucial to the design and analysis of many physical and engineering 
applications. Increases in available computational capabilities have permitted the development and testing of 
sophisticated models in the numerical simulation of turbulent flows. Direct numerical simulation, where all essential scales 
of the turbulent flow are resolved by solving the unsteady Navier-Stokes equations, are possible only at low to moderate 
Reynolds numbers. Turbulent flow analysis for engineering applications, therefore, can only be achieved by utilizing the 
time-averaged Navier-Stokes equations coupled with some level of modelling. 

The complex structure of swirling flowfields requires careful consideration of the turbulence model derivation and 
development. The analysis of turbulent transport and modeling evolves from the Reynolds -averaged Navier-Stokes 
equations and auxiliary equations for velocity and length scales for eddy viscosity specifications. Simple eddy viscosity 
models based on the Boussinesq hypothesis of linear relationship between turbulent shear stress and rate of strain have 
been quite successful in predicting a wide variety of turbulent flows. 

One of the widely used models is the two-equation k-£ model. The model developed originally by Launder and 
Spalding^ was successful in providing good predictions for a large range of turbulent flows. The equations can be derived 
from the full transport equations for the Reynolds stresses assuming fully turbulent flow. Effects such as that of rotation 
which are included in the Reynolds stress equations are cancelled out and the resulting scalar k-E equations are invariant to 
system rotation. 

For low-Reynolds number flows close to solid boundaries, adjustments to the model are needed to bridge the 
viscous dominated sublayer region with the fully turbulent flow region. The success of the wall function method depends on 
the universality of the turbulent structure near the wall. In many complex flows, however, the flow field near the wall has to 
be determined accurately and the traditional wall- function method is not satisfactory. This is because the specification of all 
turbulence quantities in terms of the friction velocity fail at separation where the flow near the wall is no longer controlled 
by the wall shear stress. Patel et al/ 6 ^ assessed the relative performance of various models which describe the near-wall 
flows and found that there are still areas of improvements needed to accurately model flow behavior near the wall. 

Jones and Launder^ extended the original k-E model to the low-Reynolds number form which allowed the calculation to 
be performed all the way to the wall. Numerical difficulties of accurately resolving the large gradients close to the wall 
necessitates resolving the wall region with very fine grid structure. Chen and Patel* 2 ) introduced a method to resolve the 
near-wall region which combines the standard k-E model with the one -equation model of Wolfshtein^ near the wall. In this 
” two-lay er" model an algebraically prescribed eddy-viscosity for the wall region is coupled to the k-E model to describe the 
details of the flow in the vicinity of the wall. Momentum and continuity equations are solved up to the wall and this reduces 
the physical uncertainties of near-wall turbulence and the numerical difficulties of resolving the vary large gradients of 
turbulence parameters. 

The purpose of this paper is to discuss the application of the k-E turbulence model with various near-wall treatments in 
the prediction of confined swirling flows. These models include, the standard wall function approach (WF), Chen and 
Patel's* 2 ) two-layer model (2L), and Lam and Bremhorst* 3 ) low-Reynolds number model (LB). 



Evaluation of the various turbulence models was performed by comparison with two selected experimental studies. The 
first is that of Daily and Nece* 4 * where rotating disk cavity circulation and secondary flows are induced by a rotating wall. 
The second is that of Roback and Johnson^ for a confined double concentric jets with a sudden expansion. Flow swirl in 
this case is induced by imposing a tangential velocity component at the outer jet. 

Numerical predictions for turbulent flows in two-dimensional axisymmetric geometries were obtained using a finite- 
volume second order upwind differencing scheme on a non-staggered grid with a pressure correction method based on the 
SIMPLE algorithm^. The development and evaluation of the turbulence models for rotating flows is part of an ongoing 
program to assess different models for rotating machinery applications. A discussion on the effects of swirl and streamline 
curvature on the turbulence structure through the gradient Richardson number formulation is given. Key problem areas will 
be identified and recommendations for the near-wall treatment as they pertain to rotating flows will be proposed. 

MODEL AND EQUATION FORMULATION 

Consider an incompressible, statistically steady and axisymmetric turbulent flow, the Reynolds averaged momentum 
and continuity equations can be expressed in a generalized form as; 

3(pd>) 3(pu<&) 1 3(pvrd>) 3 _ 3<t> 1 3 _ 

~^~ + ~^r + 7 ~ bT~ = a^ r ^a7 )+ 7a? rr ^ar )+s <t> d) 


where <& is the dependent variable 

<£>= u, v, w for the axial, radial, and tangential velocities 

p, p, and are the fluid density, viscosity and the source terms for the variable <l> 
The source terms for the dependent variables are; 


Axial direction, <b=u, T<t> x = 2p. e , T<&r ~ Pe 


5 - ^ + 


Id 


dv x 




Radial direction, <h=v, T<i> x =jJL e , npj^p^ 

„ _ v w 3p 

Sv= a7^a7 ) ' 2 ^ +p T'ar 

Tangential direction, <l>=w, r<p x = Pe, Hpr = pe 


J w“ 


pvw w 3 _ 


( 2 ) 

(3) 

(4) 


TURBULENCE MODELS In the two-equation k-£ model transport equations for the turbulent kinetic energy (k) and energy 
dissipation (e) can be written in the same general form as equation (1). 


Turbulent Kinetic energy equation 


Pt 

<I> = k, r<p x = Tct* = p + — , and S<p = G - p£ 
°k 

Energy dissipation equation 

(5) 

4 > = e, r<t>x= r < j )r = ^+ — , and S<t) = “(C 1 f 1 G-C 2 f 2 pe) 
c £ * 

(6) 


Gfc and 0g are turbulent Schmidt numbers G denotes the rate of production of the turbulent kinetic energy and is express as; 


^ ,^r,du~ dv 9 ,v- ,3u dv 9 3w- r dw w- 

C = i 2( + ^ + ( 7 )2 ] + ( - + -)2 + ( ->2, ( - _ -)2 j 


p e = p + p t , and the eddy viscosity is obtained from 


Pt “ ^p P 


£ 


Cp, Cj, C 2 , Ok; g E ^ constants whose values are 0.09, 1.44, 1.92, 1.0, 1.0, respectively. 


< 7 ) 



Near-Wall Treatment A near- wall turbulent flow can be divided into two regions, the inner viscous sublayer where low 
turbulence Reynolds number effects are important and the velocities decrease rapidly to zero at the wall, and an outer fully 
turbulent region. The successful use of the k-E turbulence model for many complex flows depends on how accurately the 
flowfield near the wall is determined. Different models are used to treat this thin sublayer region, they include: 


Wall function methods, the following equations are assumed to hold 

u + = y + fory + <11.6 

u + = ~ ln(Ey+) for y + > 1 1 .6 


where u + = ^“, y + = — — and u r = (— ) 1/>2 
^ v p 

T w is the wall shear stress which is estimated from 

= for y + < 11.6 

0-25 

* C p P u P k 5 

Tw= 5^5 for y + > 11.6 

ln(EC^ p5k°- 5 /|i) 


( 8 ) 

(9) 


( 10 ) 

( 11 ) 


where Up denotes the velocity component parallel to the wall, and 8 is the normal distance from the wall 

In this approach, k and e equations are solved with f^ = fj = f 2 = 1 only in the fully turbulent region beyond some 
distance from the wall. Boundary conditions, i.e., velocity components and turbulence parameters at that distance are 
specified in terms of the friction velocity (i^). 

In the low-Reynolds number model, the flow is resolved all the way to the wall with a very fine mesh. Many models 
have been proposed that are based on the k-E model and differ mainly in the choice of the damping functions fp, f\ and f2 to 
bridge the gap between the sublayer and the fully turbulent regions. Lam and Bremhorst’s model^ is used in this work, 
where 

V = [ 1 - exp(-0.016R y ) ]«« (1 +-^) 
fl = 1 + ( ) 3 and f2 = l-exp(-R t ) 

k 1//2 y k 2 

R v = L and R t = — are turbulent Reynolds numbers 

J V V£ 

These damping functions tend to unity with increasing distance from the wall. In order to resolve the very large gradients of 
turbulence parameters a fine mesh is required in the viscous sublayer which increases the computational time and numerical 
difficulties may be encountered. 

In order to alleviate some of the problems encountered in the low-Reynolds number approach and yet accurately resolve 
the near-wall region, Chen and Patel® pursued the two-layer concept. In this model a simple algebraically prescribed eddy- 
viscosity model for the wall region is coupled to the k-E model for the outer flow to describe the details of the flow. Unlike 
the low-Reynolds number model that requires the solution of transport equations of both k and £ all the way to the wall, the 
one -equation model requires the solution of only the turbulent kinetic energy equation in the sublayer region while 
algebraically specifying the eddy-viscosity and energy dissipation. 
v t = C(j. k I/2 L^. and e = k 3y2 /L e 

The length scales L ^ and L e contain the necessary damping effects in the near-wall region in terms of the turbulence 
Reynolds number Ry 

Lji = Ci y [ 1 - expC-Ry/A^) ] (12) 

L e = C I y[l-exp(.R y /A e )] (13) 

- 0.75 

The length scales and Lg become linear and approach Q y with increasing distance form the wall. C\ = K with Ag = 

2Q. Chen and Patel® gave values for the constant A^ = 70. The damping effects decay rapidly with distance from the wall 



independent of the magnitude of the wall shear stress. The matching between the one-equation and the standard k-e models is 
carried out along prescribed grid lines where Ry -200. 

STREAMLINE. CURVATURE AND SWIRL CORRECTIONS Turbulent flows in many engineering applications such as 
turbomachinery and combustion devices are frequently subjected to complicating influences such as mean strain and body 
forces due to rotation. In such complex flows streamline curvature and swirl can exert a large influence on the structure of 
turbulence. Bradshaw^ 10 ) reviewed the effects of streamline curvature and discussed the large effect exerted on shear-flow 
turbulence by curvature of streamlines in the plane of the main shear. So and Mellor^ 11 ^ suggested that the appropriate 

parameter governing this effect is F = , where R is the radius of streamline curvature. Militzer et al/ 12 ) provided a 

simple generalization of this parameter for a 2-D recirculating flow as 

(u 2 +v 2 ) 1/2 /R 

” (8u/3y+5v/9x) ^ ) 

They modified the turbulence production term G in the turbulent energy equations to include curvature effects and obtained 
improved predictions. Launder et al^ 13 ^ proposed a simple modification to the constant C 2 in the e-equation to account for 
streamline curvature due to swirl in the form 


C 2 = 1.0 - 0.2 R is 

where Rj s is a swirl Richardson number defined by 


w/r^ d(rw)fdT 

* s (du/9r)^ + (r 3(w/r)/3r)2 


(15) 

(16) 


Another expression 




R is = 


w 3(rw) 
Or 


(17) 


The basis of the above correction is that the effect of swirl on turbulence can be modelled through an increase in the length 
scale of the energetic turbulence eddies. 


Abujela and Lilley (14 l used a modified C 2 form (Eq. 15) with both definitions of Rj s from Equations (16) and (17) as 

applied to turbulent swirling flows. They concluded that Eq. (16) Richardson number gave better comparisons with 
experiment as compared to Eq. (17) Richardson’s number. They also found the value of C 2 obtained from Eq. (15) had to be 

limited to 0.1 ^_C 2 <.2.4 with and other constants assigned their conventional values. 


Srinivasan and Mongia^ 15 ! further split the Richardson number inlojwo parts - the swirl Richardson number R} $ and the 
curvature Richardson number R^ and corrected C 2 in the e-equations as: 


C 2 = 1.92 exp (2a s Ri s + 2a c Ri c ) 
where Rj s is given by equations (16) or (17) and 

(u 2 +v 2 ) 1/2 /R 
1 5ur d v 
r dr + dx^ 


where R is the radius of curvature given by 
0.1 and 2.4. 


R = 


(u 2 +y2)3 a 


uv( 


l3rv d u 
r dr dx^ 


(18) 

(19) 


and a s and a c are constants with values ranging between 


Chang et al^ 16 ^ investigated the streamline curvature effects in the k-e model. They managed to obtain satisfactory 
results in their hybrid k-e model where modifications of curvature effects in C 2 is made only in regions where the streamline 
curvature is large. 

In the present study curvature and swirl modifications are made to C 2 similar to Eq. (18) of Srinivasan and Mongia 
with Ri s as in Eq. (16) and R^ as in Eq. (19). The exponential form ensures that C 2 will never become negative. Numerical 

testing with several values of a s and a c reveal that C 2 may become very large and therefore, had to be limited to 0.1 <_C 2 < 
2.4. 



MODEL EVALUATION 


The various near-wall treatment models are analyzed by comparing model predictions with experimental data. Two 
cases of rotating flow experiments were selected for validation, they include; Daily and Necef 4 ) for rotating disk cavity 
experiment and Roback and Johnson^ 5 ) for swirling flow in a confined double concentric jets with a sudden expansion. The 
main criterion for selecting these cases is the different mechanisms used to generate swirling flows. In Daily and Nece 
experiment flow rotation is induced by the rotating wall, while in Roback and Johnson's Experiment, swirl is imparted to 
the flow by an outer swirling jet into a sudden expansion. Different rotation mechanisms affecting turbulence can highlight 
the differences between various turbulence models and offer certain corrections that would prove useful in accurately 
analyzing the effects of swirl. 

CASE (D - DAILY AND NECE^ In their experimental and analytic study Daily and Nece< 4 ) analyzed the steady- state 
turbulent flow in enclosed rotating disk cavities. They characterized the existence of four flow regimes depending on the 
rotational Reynolds number and cavity aspect ratio. The two-dimensional axisymmetric flow considered is that of an 
incompressible flow bounded by a disk (rotor) and a stationary end wall (stator) of a chamber as shown in Figure 1 . The ratio 
of the axial clearance between the rotor and the stator (s) to the radius of the disk (a) is 0.0255. The disk rotates with a 
rotational Reynolds number R=4.4xl0^ defined as R=fta^/v, where ft is the disk rotational speed in rad/sec and v is the 
kinematic viscosity. 

Numerical computations were performed on a 33x75 grid with different grid clustering near the walls for the different 
near- wall models. Figure 2, shows the velocity vectors at the top region of the cavity using the WF model. Centrifugal 
forces move the fluid radially outward on the disk, axially away from the disk on the wall casing, and radially inwards on the 
stationary end wall. Figure 3, shows the axial variations of the radial velocity component (v) at a radial position r/a=0.765. 
The agreement is fair with some discrepancy for all near-wall models close to the rotating disk. Figure 4, shows the axial 
variation of the tangential velocity (w) component at the radial position. At the rotating disk (x=0), the tangential velocity 
component approach the value (aft). The 2L near- wall model seem to offer closer agreement with the data than the other two 
models. 

The presence of comer regions presents a difficulty in defining the normal distances used in the definition of turbulent 
Reynolds number (Ry). In the present analysis, values of the normal distance from a wall were based on the normal distance 
to the nearest solid boundary. Streamline curvature and swirl corrections have not been used in this case. 

CASE (21- RQEACK AND JOHNSON ^ Predictions of the experiments of Roback and Johnson^ have been presented by 
several workers, e.g. Sloan et a]/ 17 ) and Durst and Wennergerg* 18 ). Unfortunately, inlet profiles were not provided in their 
experiment. Therefore, calculations were started at the expansion plane using the measured velocity profiles at 5mm 
downstream of the expansion after some adjustments near the edges of the coaxial jets. Measurements of all three main 
turbulent intensities were used to calculate inlet values of the turbulent kinetic energy. Energy dissipation rate was 
estimated from 

Cuk^ 

e =— ^ — (20) 

where L is a length scale of turbulence at the inlet of the order of L=10"^ m. 

Figure 5, shows an illustration of the test chamber geometry. The confluence plane of the primary (inner) and 
secondary (outer) jet streams coincides with the chamber expansion plane. The chamber diameter is about twice the 
secondary tube diameter. The exit from the 8-bladed, 30°, free vortex swirl generator is located approximately 0.05 m 
upstream from the confluence plane. 

A prevalent phenomenon in axisymmetric swirling flows in such geometries is the "bubble" or vortex breakdown 
which has been studied extensively* 19 * 20 ’ 21 ’ 22 ). The near axisymmetric breakdown can be partially understood from a 
simplified analysis of the role of pressure and centrifugal forces. It is identified by a slowly varying vortex core which 
undergoes an abrupt and rapid deceleration, forming a free stagnation point, followed by a region of flow reversal. It is 
known that the structure of vortex breakdown is unstable and asymmetric in the azimuthal direction, and displays 
unsteadiness in the axial direction* 23 * 24 ). However, no periodic or nonaxisymmetric behavior attributable to the vortex 
breakdown was observed in Roback and Johnson's experiment. 

In the numerical simulation of the experiment, a 150x100 grid nodes was used with different clustering on the walls for 
the different near-wall models used. Figure 6, shows the velocity vectors indicating the presence of a closed recirculation 
zone at the center with additional zones at the comer downstream and between the inner jet and the outward diverted 
secondary jet. The figure also shows flow diversion outwards with high gradients characterized by large turbulent shear and 
fluctuation levels. 

Comparisons were made of the radial variations of flow variables at two axial locations, x=0.025 m upstream of the 
vortex bubble and x =0.1 02 m located inside the vortex bubble. Figure 7a, shows the radial variation of the axial velocity 



profile at x=0.025 m using the WF method, 2L model and LB model. Fair agreement is predicted by the different models. 
They also seem to predict small negative velocities at a radial position r-0.0153 m (the interface between the inner and 
outer jets), slightly underpredicting its strength and width. Figure 7b, shows the radial variation of the axial velocity 
profiles at x=0.102 m. The 2L model shows a better agreement with the experimental data. These velocities are slightly 
underpredicted above the outer jet diameter. 

Radial variations of the tangential velocities are shown in Figure 8a and 8b at x=0.025 m and x=0.102 m respectively. 
Figure 8a shows that the 2L model offers a better agreement with the experiment as compared with the WF and LB models. 
At x=0.102 m. Figure 8b shows that the swirl velocity is underpredicted. That is because the radial transfer of 
circumferential velocity is highly dependent on the turbulent diffusion mechanisms which are not accurately modelled in the 
isotropic eddy-viscosity k-e model used here. 

The turbulent intensity predictions for the k-e model using the different near-wall treatments seem to follow similar 
trends as shown in Figures 9(a,b), 10(a,b), and, ll(a,b). In general within the approximations of the isotropic k-E model, 
the 2L model offer a marginal improvements over the WF and LB near wall models. The peaks in the axial, radial, and, 
tangential turbulence intensities occur around the edges of the inner and outer jets. Figure 13a, shows the axial -azimuthal 
Reynolds stress profile at x=0.025 m. Figure 12b and 12c, show the axial -radial and radial -azimuthal Reynolds stress 
profiles at x=0.102 m. 

The analysis of the main turbulent intensities and of the Reynolds stress components using the isotropic eddy- 
viscosity k-£ turbulent model do not reveal exclusively the advantage of one near-wall model over the other. Moreover, 
Reynolds shear stress profiles are sensitive to the upstream inlet conditions and the developing mean flowfield. Although 
the mean flow quantities show a general trend of improved predictions using the 2L near-wall model, the main effects of 
turbulence are due to anisotropy of Reynolds stresses especially around the highly sheared region of the outward diverted 
outer jet and the vortex bubble. 

Streamline curvature and swirl corrections have been attempted in the present analysis with little success. Corrections 
of C 2 using equation (18) with equations (16) and (19) for the swirl and curvature Richardson numbers. Figure 13(a,b) show 
the radial distribution of the radial velocity at x=0.025m, and the axial velocity at x=0.102m. Small improvement is 
detected with these corrections. The constants a, and used are those recommended by Srinivasan and Mongia^ 15 ^ in their 
calculations (ct s =-0.75 and a c =-2.0). These constants were not optimized in the present calculations. 


CONCLUSIONS 

Flow predictions were performed for the standard k-e turbulence model with different near-wall treatments to assess 
their performance when applied to rotating flows. Comparisons of predictions with the experimental data of Daily & Nece, 
and Roback & Johnson show reasonable agreement for all near-wall models and in general, the two-layer model seem to 
offer better comparisons compared to the wall function and Lam & Bremhorst low -Reynolds number models. From a 
computational perspective, the two-layer model require less computer time and relatively few grid points in the wall region 
than the low -Reynolds number model and is less sensitive to the location of the interface between the sublayer and the fully 
turbulent regions. Streamline curvature and swirl corrections show small improvements. However, further study is needed to 
optimize their constants. 
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Fig. 9a Radial Distribution of Axial Turbulence 
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ADVANCED TURBULENCE MODELS FOR TURBOMACHINERY 


Ali H. Hadid, Michele E. DeCroix, and Munir M. Sindir 
Rocketdyne Division, Rockwell International 


ABSTRACT 

Development and assessment of the single-time-scale k-e turbulence model 
with different near-wall treatments and the multi-scale k-e turbulence model for 
rotating flows are presented. These turbulence models are coded as self- 
contained module decks that can be interfaced with a number of CFD main flow 
solvers. For each model, a stand-alone module deck with its own formulation, 
discretization scheme, solver and boundary condition implementations is 
presented. These satellite decks will take as input (from a main flow solver) the 
velocity field, grid, boundary condition specifications and will deliver turbulent 
quantities as output. These modules were tested as a separate entities and 
although many logical and programming problems were overcome only wider 
use and further testing can render the modules sufficiently ’’fool proof". 
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STRUCTURE OF FUTURE CODES LEADS TO 
DEVELOPMENT OF MODULES 
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TURBULENCE MODEL DECK STRUCTURE AND 
INTEGRATION WITH NAVIER-STOKES SOLVER 
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SINGLE-SCALE k-e MODEL (CONT’D) 
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STRONGLY IMPLICIT SOLVER 


KEMOD-1 MODULE DECK (CONT'D) 
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THE SUBROUTINE SOLVES THE DISCRETIZED EQUATIONS AFTER 
MODIFYING THE SOURCES AND BOUNDARY CONDITIONS FOR THE 
PARTICULAR PROBLEM 
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MOMENTUM EQUATION SOLVER OF THE MAIN ROUTINE. IT UPDATES 
THE FLUX SOURCE TERM OF THE DISCRETIZED MOMENTUM EQUATION 
DUE TO WALL SHEAR STRESSES 
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KEMOD-1 MODULE DECK (CONT’D) 
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MULTI-SCALE k-£ MODEL (CONT'D) 
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SUBROUTINE CALCE ASSEMBLES THE COEFFICIENTS AND SOURCE 
TERMS FOR THE DISCRETIZED K p , e p , 1^, e t TRANSPORT EQUATIONS 
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Abstract 

A finite-volume procedure is used to 
compare the performance of different high order 
turbulence models for confined swirling turbulent 
flows. Eddy-viscosity single and multi-scale k-e 
turbulence models together with second r moment 
algebraic and Reynolds stress closure models are 
tested for a two-dimensional, axisymmetric swirling 
flow case. The ability of second-moment closure 
models to capture the interaction between swiri and 
the turbulent stress field is crucial to the predictive 
performance of the computational scheme. 

To enhance the predictive capability of CFD 
tools for engineering applications, advanced 
turbulence models are coded as self-contained module 
decks that can be interfaced with a number of CFD 
solvers. Three of these modules, namely the single 
and the multi-scale models and the Algebraic stress 
model (ASM) have been successfully interfaced and 
tested with the code MAST of the University of 
Alabama at Huntsville in a relatively short time. 
These modules are independently tested and evaluated 
with the data of Roback and Johnson for swirling 
turbulent flow in a confined double concentric jets 
with a sudden expansion. 

Modularization of a general purpose CFD 
code structure in terms of different aspects of physical 
models is necessary for computational efficiency. 
Further, individual modular routines are transportable 
and can be easily modified to include extra physical 
effects. This would allow many users using different 
CFD codes to concentrate their talents on developing 
and improving physical hypothesis for specific 
engineering problems. 

Introduction 

Computational Fluid Dynamics (CFD) has 
been used extensively for the last decade or so in 
analyzing complex flow phenomena for many 
industrial applications, such as, turbomachinery and 


combustion devices. Most flows of technological 
interest are turbulent and for many of them, relatively 
simple prediction methods are sufficient to produce 
results of engineering accuracy. For others, mainly in 
high technology applications, accurate predictions 
using high order turbulence models are required. 
Increases in available computational capabilities have 
permitted the development and testing of sophisticated 
models in the numerical simulation of turbulent 
flows. Direct numerical simulation, where all 
essential scales of the turbulent flow are resolved by 
solving the unsteady Navier-Stokes equations, are 
possible only at low to moderate Reynolds numbers. 
Turbulent flow analysis for engineering applications, 
therefore, can only be achieved by utilizing the time- 
averaged Navier-Stokes equations coupled with some 
level of modelling. The analysis of turbulent 
transport and modelling evolves from the Reynolds- 
averaged Navier-Stokes equations and auxiliary 
equations for velocity and length scales for eddy 
viscosity specifications towards a more sophisticated 
modeling startegy - one offering greater width of 
applicability, particularly in complex shear flows or 
where external force fields modify the turbulence 
structure. 

One of the widely used models is the two- 
equation single- time-scale k-e model®. In this model 
transport equations for the turbulence energy (k) and 
the energy dissipation (e) are solved to determine the 
turbulent eddy viscosity. An improvement to the 
single scale k-e model is the multi-time-scale k-e 
model where the energy spectrum of a turbulent flow 
is split into a production range and a dissipation 
range®. Improved predictions using the multi-scale 
over the single-scale k-e model have been 
demonstrated® 4 ^. 

Other complicated single-scale models 
offering greater width of applicability, particularly in 
complex shear flows or where external force fields 
modify the turbulence structure are based on secon- 
moment closures. These take the exact equations for 

the transport of the Reynolds stresses (ujuj ) as then- 
starting point and devise approximations for the 
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unknown turbulent correlations appearing in them. In 
a three-dimensional flow, or even in an axisymmetric 
flow, all six components of the Reynolds stress 
tensor are nonzero. With a full second-moment 
closure model (RSM), therefore, differential transport 
equations need to be solved over the solution domain 
for each of these components. This represents an 
increase in the task of numerical solution compared 
with the situation where the k-e eddy-viscosity model 
is adopted. An intermediate level of modeling has 
evolved^* 6 ) known as Algebraic second-moment 
closure (ASM), with the aim of retaining the greater 
physical realism of second-moment treatments while 
achieving computational times closer to that of an 
eddy-visocity model. The simplification is achieved 
by approximating the convective and diffusive 
transport of the Reynolds stresses in terms of the 
corresponding transport of turbulent energy. This 
allows the transport equations for the stresses to be 
expressed as a set of algebraic formulae containing the 
turbulence energy and its rate of dissipation as 
unknowns. Second moment schemes have been 
extensively and successfully applied to a wide range 
of flows, as reviewed for example by Leschiziner^ 7 ). 
Few applications, however, have considered 
axisymmetric swirling flows^ 6 * 8 ) where the external 
forces due to swirl exert damping effects on the 
turbulent transport 

Progress in turbulence modelling have been 
paralleled by improvements in numerical techniques, 
essentailly, combining second moment closure with 
non-orthogonal, co-Iocated grids using finite-volume 
methods. However, the implementation of RSM into 
non-orthogonal finite-volume codes poses difficulties: 
the co-located variable arrangement can cause 
decoupling of the mean velocity and Reynolds stress 
fields leading to oscillating solutions or even 
divergence. Using a special interpolation procedure in 
the context of Rhie^\ Obi and Peric^) calculated 
the two-dimensional turbulent flow on a co-located 
grid arrangement usin the Reynolds stress turbulence 
model. 

In the present paper, we pesent predictions of 
two dimensionaVaxisymmetric swirling flow using 
various models based on eddy-viscosity single and 
multi-scale k-e and on second moment closure. These 
models are cast in a modular form enabling them to 
be used with a number of flow solvers based on the 
finite-volume and finite-difference methods. A 
discussion of the different models used and their 
assessment is presented. The modular structure of the 
different turbulence models will also be presented and 
discussed. 


Theory and Model Equations 


Reynolds averaged continuity and momentum 
equations which may, respectively, be written as 


MZ + i2erV = 0 

dx r 9r 


(1) 


aoruo aprvo a , a<t\ a , a<t>, „ 
“ + ar _= ^ (r ^ar >+ * (r ^ar )+rS<I) (2) 


Where <t> stands for any of the momentum 
components U, V, and rW and the corresponding 
sources S<j> are 


Su = 


dP dpu 2 1 3pruv 
3x dx r 3r 


„ _ _ 3P rW 2 2uV J_ 3prv 2 3puv nv- 
V ~ 3r + r r 2 r Jr 3x + r 


2u 3rW 3(puw) 3fpvw) 


S ' w= -Tir' r *r- r 


dr 


- 2pvw 


where p,p. are the fluid density and visocosity 
respectively. 

The appearance of the Reynolds stresses tijui 
represents an unknown correlation and different 
turbulence models provide the means of relating these 
unknowns to known determinable quantities. 

Single-Scale Eddv-Viscositv Turbulence 
Models 


Here it is assumed that a single-time-scale 
(proportional to k/e) can be used to describe the 
turbulent flow. Turbulence is simulated through 
transport equations for the turbulent kinetic energy 
(k), and its rate of dissipation (e). The stress tensor is 
modelled using a gradient transport model of the form 
- 3Ui 3U; 2 

- u iuj = v t(^7+3^)-jk8 ij « 

The generalized form of the two-equation eddy- 
viscosity turbulence model can be written as 
Kinetic Energy (k) equation: 


wnere 


^ UKJ I 

Ck_ 3xi 


Ck = Dk + P - e 

aUjk 


_ a r . v t 3k. 
D k=^t( v +-)^] 


ok'ax; 


(4) 


Convection of k 


Diffusion of k 


The turbulent flow considered is two- 
dimensional and steady which can be described by the 
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k 2 

V[ = eddy viscosity = Cjj. — 


Energy Dissipation (e) equation: 


C e = D e + £(CeiP-C e2 e) 


(5) 


where 


C e = 


9U;e 


9x 


J 


T-, 9 r , Vt 9t 


Convection of z 
Diffusion of z 


In the present study, the standard two-equation model 
was used with the wall function^ 1 ) and the two-layer 
model(H) to bridge the gap between the near-wall 
log-layer region and the fully turbulent region away 
from the wall In the standard model the numerical 
values of the constants are C^O.09, C e i = 1.44, 
C e 2= : l-92, Oit=I. 0 and cr e =1.3. Details of the 
implementation of the wall function and the two-layer 
models can be found in Hadid and Sindi/ 12 ). 

Multi-Time-Scale k-e Turbulence Model 

The Multi-time-scale turbulence model used 
here is based on the variable energy partitioning of 
the turbulent energy spectrum proposed by Kim and 
Chen^). In this model the turbulent kinetic energy 
spectrum is divided into two sets of wave number 
regions giving two evolution equations for each 
region. These equations represent the kinetic energy 
(k p ) and the energy dissipation (e p ) in the production 
range of the spectum and the kinetic energy (kt) and 
the energy dissipation (eO in the dissipation range of 
the spectrum. This model allows the partition to 
move toward the high wave number region when 
production is high and toward the low wave number 
region when production vanishes. 

The equations for the turbulent kinetic energy (k p ) 
and the energy transfer rate (e p ) for the production 
range are 


Ckn - E>k D + P - £r 


e. 


Qt = D et + t ( c tl £ p + Ct2 e t) - Ct3 r- 
K t M 


■ _ 3Djk 0 

where Cfcp = p ^ / and 


( 9 ) 



K P - axj LV>i aiq/flXj 5 


Dkn = and 




similarly for Cfc t , c £t» and D £t equations and 
the model constants used are those of Kim and 
Chent 3 )* 

2 

i p 2 £ p 

The terms (- C p i — ) and (p C 1 1 — ) 

p ^ Kp k t 

represent variable energy transfer functions. The 
former increases the energy transfer rate when 
production is high and the latter increases the 
dissipation rate when the energy transfer rate is high. 
The turbulent viscosity mt is given as pt = 

pCjifk^/Ep = pCjjk*Vet, where k=k p -rk t is the total 
turbulent kinetic energy. 

Second Moment Gosure Models 

The exact form of the Reynolds stress 
equations can be derived from the time-averaged form 
of the Navier-Stokes equations and can be written as: 


D 


£ t (pujuj) = Py -M^ij + Dy +pey 

where 

aui _ aui 

Pij = -p c + uiufc fei } 


(10) 


Production 


( 6 ) 


D'j = 9xk ( P u ‘ u j uk ‘ ^ pui ' pu j ) Diffusion 

<p;; = p Pressure-strain redistribution 

-> oxj OXj 


C ep - D ep + p k p ^pl p + C p2 £ p) - C p 3 jT (?) 

The equations for the turbulent kinetic energy (kt) and 
the dissipation rate (eO for the high wave number 
transfer region are 


e ij - 3 §ij e 


Dissipation 


Due to the introduction of correlations of higher 
orders, modelling of these terms is required to close 
the set of equations. 


Ckt - Dkt + £p - £t 


( 8 ) 
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Algebraic Stress Model (ASM) 


The ASM model used is based on the work 
of Rodtf^. The idea is to simplify the stress equation 
(eq. 10) by approximating the convective and diffusve 

transport of the Reynolds stresses (upj) in terms of 
the corresponding transport of turbulent energy. This 
simplification allows the transport equation of the 
stresses to be expressed as a set of algebraic formulae 
containing the turbulent energy and its rate of 
dissipation as unknowns. This set of algebraic 
equations can be written as; 

u « u j = f Pi j ' f 5 »j e ^ij 1 (11) 

The pressure-strain term is decomposed into a 
fluctuating part (<by t i), a part due to the mean rate of 
strain and a part due to reflected wall- 

influence (<*>ij iW ), i-e., <t > ij = <Dy # i + + <&ij,w 

Rotta’s return to isotropy concept is used to 
model the non-linear pan (Oyj) as 

^ij.l = - C lf ( u i u j - fkSij ) 

<X>ij f 2 is modelled using the isotopization of 
production concept as 

<&ij,2 = - C2 ( Pij - |p5jj ) 

The wall reflection term <t>y jW is modelled 
following Shir^ 1 ^ and Gibson and Launder^ 14 ) as 
^ij,w = ^ij.lw + ®ij f 2w 

where 
^ij t lw = 

* £ — 3 — 3 — 

p J(uk^m n k n m 5 ij -J^i n k n j “ 2 u k u j n k n i ) f (12) 
^ij,2w ~ 

* 3 3 

^2 (^km.2 n k n m^ij"2 ^ik.2 n k n j _ 2 ^jk.2 n k n i) f (13) 
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B = ^ + 2(1^2) (^H’lw + ^1 l,2w) 
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p 
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where n* is the wall-normal unit vector in the i- 
direction. The wall distance function (f) represents the 
ratio of turbulence length scale and the wall distance 
C« 75 k l-5 , 

m 1 

f = ( ) — where An is the wall-normal 

K £ An 

distance. The above wall-correction terms are written 
in a tensorialy invariant form and their effect is to 
transfer energy from the wall-normal normal stress 
component to the tangential stresses i.e it is 
redistributive. 


Reynolds Stress Model (RSM) 


In the RSM model the full transport equation for the 
Reynolds stresses (eq.10) are solved for each stress 

component (u[uj) after modelling the diffusion and the 
pressure strain terms similar to Launder et. 

The diffusion term is modelled as 


_ 3 r _ — k auku] 


For ax i symmetric swirling flows the set of 
algebraic stress equations can be written in a general 
matrix form as A T = B where 


The pressure-strain redistribution term is 
modelled in a similar way to that used in the ASM 
model discussed earlier. Special consideration is given 
to the problem of mean velocity-Reynolds stress 
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decoupling which appear when using a collocated grid 
arrangement which is a source of numerical 
instability. This is done by invoking a special 
interpolation procedure for the cell-face stresses in the 
context of Rhie(9). This practice results in the 
addition of normal stresses to the pressure term where 
the cell- face velocity is sensitized to the pressure 
differences as well as to normal stress diffrences at the 
nodes surrounding the face. 


Turbulence Model Decks (Modules) 

As the state-of-the-art of computers has 
advanced, so has the range, size and complexity of 
flow models being applied. Users have become more 
sophisticated and there is a constant demand for 
improvement. CFD codes have adapted to this 
demand and many general-purpose computer codes 
have been developed and used. As general-purpose 
codes become larger, their code structure becomes 
sophisticated. In general codes can be divided into 
three main areas, they include; 1) Numerical 
algorithms (which can be subdivided into 
discretization methods and solution techniques). 2) 
Methods of dealing with complex geometries. 3) 
Physical models (which include turbulence models, 
porosity, combustion kinetics, two-phase flow...). It 
seems, therefore, that the practicing engineer must 
have the knowledge of all these elements of the CFD 
program in order to successfully utilize this code. To 
obtain the maximum benefits from these general- 
purpose CFD codes, modularization of the code 
structure may be necessary. That is developing 
individual modular routines for the solver and for 
different physical models for example. If such 
modules are successful it would allow users to 
concentrate their talents on developing and improving 
physical hypothesis such as turbulent models for 
example that can easily be tested using such modules. 

In the present work, turbulent modules are 
being developed to meet this need. Figure 1, shows a 
flowchart of a turbulence module interfaced with a 
typical main flow solver. The module is called by the 
flow solver passing to it the mean flow velocities, 
mass fluxes at cell faces and grid information among 
others. The turbulence differentail equations are 
discretized and the matrix coefficients are setup and 
solved using Stones strongly implicit method^ 16 ). In 
the ASM module, the set of algebraic stress equations 
are solved simultanously using Gauss-Seidel method 
at each step or iteration. In the eddy-viscosity models 
the values of k, e, and eddy viscosity (p.[) are passed 
to the main flow solver, while, in the second moment 

closure models the Reynolds stresses lijuj are passed 
to the main solver. The solver then calls subroutine 
MODIFY of the module where the momentum 


sources are modified to account for the near- wall shear 
stresses in the eddy-viscosity models or to calculate 
Reynolds stress gradients in the second moment 
models. 

These modules are structured to be self- 
contained and transportable to a number of general 
purpose CFD solvers to maximize computational 
efficiency. They have been tested independently at the 
University of Alabama at Huntsville using the 
MAST code. 

Results 

The various turbulence models are analyzed 
by comparing model predictions with the 
experimental data of Roback and Johnson^ 1 ?) for 
swirling flow in confined double concentric jets with 
a sudden expansion. 

Figure 2, shows the decay of the mean axial 
centerline velocity using both the single and multi- 
scale k-e models. Figure 2a, shows the comparison 
using the wall-function near wall approach and Figure 
2b shows the results using the two-layer near wall 
model. The single-scale k-e model seems to 
underpredict the extent of the central recirculation 
zone as compared with the multi-scale k-e model. 
Moreover, improved comparisons with the data are 
obtained using the two-layer near wall model. Figure 
3, shows the radial profile of the mean axial velocity 
at two distances downstream of the jet exit. Again, 
the two-layer model predicts better comparisons with 
tha data than the wall function approach. The radial 
profiles of the mean tangential velocity are shown in 
Figure 4. 

Figure 5, shows the radial profile of the 
mean axial velocity at three axial locations using the 
algebraic stress model (ASM) with wall function and 
two-layer near wall models. The radial profiles of the 
rms axial turbulent intensity are shown in Figure 6. 
Streamline contours are shown in Figure 7 using the 
single-scale k-e and the ASM models with the two- 
layer near wall approach. The extent of the central 
vortex is better predicted using the ASM model. 
Preliminary results were also obtained using the full 
Reynolds stress model (RSM). Comparisons with the 
backward facing step data of Driver and 
Seegmiller(^) shows improved predictions over the 
single scale k-e model as shown in Figure 8 where 
the radial profiles of the axial normal stress and shear 
stress are plotted at four step heights downstream. 
Further testing of the RSM model for swirling flows 
are planned. 
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Conclusions 


Different turbulence models for industrial 
applications have been formatted in a modular form 
and successfully interfaced and tested independendy 
using two different main flow solvers. The turbulence 
models include the single and multi-sclae k-e models 
both with wall functions and two-layer near wall 
models. Second moment models that include the 
algebraic (ASM) and full Reynolds stress model 
(RSM) have been tested. It was shown that the two- 
layer near wall model improves predictions as 
compared to the wall function approach. Convergence 
of the stiff ASM model equations was obtained by 
solving the 6x6 stress equations (for 
axisymmetric/swirling flows) at each iteration. The 
wall-reflection terms in the pressure-strain model 
showed little or no improvements in the ASM model 
predictions. Elaborate pressure-strain models that 
require no wall-damping are needed e.g. Speziale 
et.al^). The full Reynolds stress model (RSM) 
promises to be the next model to be used for 
engineering applications. 
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Figure 2. Mean axial centerline velocity 
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Figure 6. Radial profiles of the axial turbulent intensity 
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Figure 7. Sreamline contours 




Figure 8. Backward facing step (Driver & Seegmiller) 
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A NUMERICAL STUDY OF TWO-DIMENSIONAL VORTEX 
SHEDDING FROM RECTANGULAR CYLINDERS 

A. H. HADID * M. M. SINDIR* R. I. ISSA 1 

Abstract 

An efficient time-marching, noniterative calculation method is used to analyze time- 
dependent flows around rectangular cylinders. The turbulent flow in the wake region 
of a square section cylinder is analyzed using an anisotropic k-£ model. Initiation 
and subsequent development of the vortex shedding phenomenon is naturally cap- 
tured once a perturbation is introduced in the flow’. Transient calculations using 
standard eddy-viscosity and anisotropic k-c models, averaged over an integral num- 
ber of cycles to get the fluctuating energy (organized and turbulent), are compared 
with experimental data. It is shown that the anisotropic k-c model resolves the 
anisotropy of the Reynolds stresses and gives mean energy distribution closer to the 
experiment than the standard k-c model. 

1. INTRODUCTION 

Vortex shedding is a periodic unsteady flow phenomenon that occurs frequently behind bluff 
bodies and is therefore of great practical importance. Many attempts to calculate the two- 
dimensional (2-D) vortex shedding motion past square and .circular cylinders by solving the 
unsteady Navier-Stokes equations w^ere successful at low Reynolds numbers where the flow is 
laminar and the fluctuations are periodic, e.g., [l] and [2]. At higher Reynolds numbers which 
are more relevant in practice, turbulent fluctuations are superimposed on the periodic unsteady 
motion. The problem then concerns the decomposition of the flow into organized motion that is 
resolved in the calculation and a remaining turbulent motion to be represented by a turbulence 
model. Previous analysis of vortex shedding calculations at high Reynolds numbers have not 
been successful due to the inadequacy of the standard k-c model and the lack of affordable 
higher order models that take into account the anisotropy of the turbulent intensities 

Franke et ai [3] analyzed the unsteady turbulent flow for a square cylinder using the stan- 
dard k-c model. They showed that the model tends to damp the periodic shedding motion 
underpredicting the Strouhal number. They also analyzed the detailed experimental results of 
Cantwell and Coles [4] for vortex shedding in the 2-D wake behind a circular cylinder. They 
additionally point out the need for improved models that account for the history and transport 
effects of the individual stresses. Maclnnes et ai [5] used the standard k-c model to simu- 
late the periodically forced turbulent mixing layer investigated experimentally by Weisbrot and 
Wygnanski [6]. They managed to capture the main features of the mixing layer development 
where there is a clear distinction between the organized and the random turbulent motion. 
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Majumdar and Rodi [7] have shown that the separated turbulent flow past circular cylinders 
cannot be predicted realistically without a time-accurate numerical procedure to account for 
the periodic shedding of vortices. 

Experimental investigations are needed to judge the different numerical and turbulent 
schemes. Durao et ai [8] conducted an experimental study of transient turbulent flow be- 
hind a square cylinder. They used spectral analysis and digital filtering of the LDV data in 
order to separate and quantify the turbulent and periodic, nonturbulent motions. They show 
for example that in the zone of highest velocity fluctuations the energy associated with the tur- 
bulent fluctuations is about 40 % of the total energy. Therefore, for a successful simulation of 
transient turbulent flows, a reliable time-accurate numerical procedure and a good turbulence 
model are needed. 

The purpose of the present paper is to model turbulent vortex shedding flows using an ef- 
ficient time-accurate numerical procedure based on the PISO [9] methodology. Calculations of 
the turbulent vortex shedding are performed using the two-equation k-e model with isotropic 
eddy- viscosity and with a modified two-equation model using an anisotropic eddy-viscosity. In 
the anisotropic model, nonlinear corrections are added to improve the eddy-viscosity represen- 
tation of the Reynolds stresses as developed by Yoshizawa [10] with the aid of a two-scale direct 
interaction approximation. A similar anisotropic eddy-viscosity model was also developed by 
Speziale [11]. The adequacy of the models to simulate transient turbulent flows is assessed 
with the aid of the experimented results of Durao et al. [8] for vortex shedding in the 2-D wake 
behind a square cylinder at Re = 14,000. 


2. MODEL EQUATIONS 

The basic equations of motion in transient periodic flows can be written after separating the 
flow into an organized (phase averaged) component 

1 N 

Ui(xi,t) - — i + nT) (1) 

iV n-0 


where U t (x t) i) is the resolvable portion of the instantaneous velocity u t , and T is the period 
of the oscillation, and a random turbulent component v! x (x u t). The instantaneous velocity 
u t (xi,t) is then given by 

u t = U % -f u\ = u x 4- u'i (2) 


where is the time-mean component of the velocity, and TT X is the periodic fluctuating compo- 
nent. Assuming an incompressible flow, the momentum equations can be written after applying 
phase averaging as; 


8U t 


+ U, 


dU, 

d Xj 


1 dP 

p dx t + 


d dUj 

d Xj dxj 
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( 3 ) 


where R XJ = — < id, u* } > is the phase- aver aged Reynolds stress tensor and v is the kinematic 
viscosity. 

Standard Isotropic k-e Model 

In the standard isotropic k-e model [12], R tJ is approximated by using the eddy- viscosity i/ t as; 
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where k is the phase- aver aged turbulent kinetic energy and v t = C^(k 2 /e), e is the phase- 
averaged energy dissipation rate and C ’ M is a model constant. The spatial and temporal distri- 
bution of k and e are determined from differentia] transport equations of these quantities 
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( 6 ) 


where G — R tJ *— — * is the turbulent generation term. The constants C fll Cj, C 2f 

d Xj 

have values of 0.09, 1.44, 1.92, 1.0, and 1.3, respectively. 

Anisotropic k-e Model 

In the anisotropic model the Reynolds stresses can be expressed as; 


cTfc, and <7 e 


R,j = —~kS tJ A v t 
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( 11 ) 


and Cr m (m = 1, 2, 3) are model constants. The first two terms on the right hand side of 
(7) give the familiar isotropic eddy- viscosity representation, while the third and fourth terms 
express the anisotropy of These additional nonlinear quadratic terms of the mean velocity 
gradients seem to be a simple way to resolve the individual normal stresses with the k -e model. 
The anisotropy is reflected especially in the k -e equation where both the diffusion and pro- 
duction terms are quadratic forms of the mean velocity gradients and turbulent kinetic energy 
gradients. 

The anisotropic eddy- viscosity model has been successfully used by Nisizima and Yoshizawa [13] 
and Myong and Kasagi [14] for fully developed turbulent channel flows. In their calculations 
only Cti and Cr 2 were optimized to reproduce the anisotropy of the turbulent intensities since 
Ct 3 does not appear in their equations. In the present study the flow is shear dominated with 
little departure from isotropy. Therefore, the model constants Ctj, CV 2 , and Cr 3 were opti- 
mized to 0.01, 0.01, and 0.001, respectively, to satisfy the realizability constraint. (Note: zero 
constants reduce to the isotropic eddy-viscosity model.) 

Applications of the k -e isotropic and anisotropic eddy- viscosity models were made using 
wall functions to bridge the viscosity affected near the obstacle wall region. It is assumed that 
inadequacies in near-wall modelling play a minor role to the inaccuracy of normal Reynolds 
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(a) Normal velocity fluctuation along the (b) Power spectrum of the normal velocity 
centerline fluctuation 



(c) Streak line plot 

Fig.l: Flow characteristics of the wake for Re = 14,000 


stress differences arising from use of an isotropic eddy viscosity. Improvements can be made by 
integrating all the way to the wall [15] or by using the two-layer model of Chen and Patel [16]. 


3. NUMERICAL METHOD 

The PISO methodology [9], in conjunction with a finite-volume technique, is used to solve the 
implicitly discretized, time-dependent flow equations. The method is essentially noniterative, 
where the solution process is split into a series of steps whereby operations on pressure are 
decoupled from those on velocity at each time-step. The avoidance of iterations substantially 
reduces the computational effort compared with that required by iterative methods. This is 
possible since the splitting error of PISO is negligibly small at the level of time-step required 
to eliminate the temporal truncation error. A backward temporal difference scheme is used, 
while the convective terms are discretized using a second-order upwind difference scheme. The 
method can also be used for steady-state flows, e.g., Hadid et ai [17]. 

Calculations are performed for the turbulent flow around a square cylinder (step height, 
H = 20 mm) in a domain extending about 16 H downstream and 2.5 H upstream of the obstacle. 
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Fig. 2: Centerline distribution of mean axial Fig.3: Time-mean kinetic energy of the ve- 
velocity locity fluctuations 

The calculations captured the vortex shedding phenomenon after perturbing the flow at the 
inlet. A reference velocity of 0.68 m/s and turbulence intensity of 6% (i.e., k =< u' 2 > = 
3.6 x 10” 3 m 2 /s 2 ) were used as the inlet conditions. The length scale L of turbulence at the inlet 
was not measured in the experiment but an order of L ~ 0.1 mm was assumed from which the 
energy dissipation rate e = k z t 2 /L was estimated. It is expected that the calculated results are 
not sensitive to the precise value of e used at the inlet. The upper and lower boundaries are 
treated as symmetry planes, at the exit, a zero-gradient outflow boundary condition is applied 
to each variable. The computational domain is resolved by 75x40 grid cells with clustering at 
the obstacle walls. An optimized time step of 0.001 sec. was chosen for the calculations. 

4. RESULTS AND DISCUSSION 

Figure 1(a) shows the normal velocity history at the centerline of the wake for Re = 14000 at five 
step heights downstream. The power spectrum of the normal velocity fluctuations (Fig. 1(b)) 
confirms the oscillatory nature of the flow with a single predominant frequency of about 4.7 Hz, 
which is in agreement with experimental results [8]. Figure 1(c) shows a marker particle trace 
at time=3sec., which illustrates the shedding pattern. In order to calculate the time-mean 
kinetic energy of the velocity fluctuations, the fluctuating velocity component (organized -f 
turbulent) is u, = u x — U t . For the 2-D plane geometry considered, the kinetic energy of the 
velocity fluctuations can be written as; 

E ~ 4 (^i 2 + “a*) ( 12 ) 

where uC = u 2 — 2 u l U t + U t 2 } and the time-mean value of the kinetic energy of the velocity 
fluctuations is 


3 
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o Dura o et al. [8] 

— Anisotropic K-€ Model 

— Anisotropic K-€ Model with Zero 
Production of K Upstream of Obstacle 




Fig. 4: Centerline distribution of mean axial Fig.5: Time-mean kinetic energy of the ve- 
velocity locity fluctuations 

where 0 

s ; 2 - (u 2 - 2 ujTi + it: 2 ) = [(Ui + - 2 (u x + u\) Ui + u:\ 

and from the definition of time averaging u[U x = 0, we get, 

(*=1,2) (14) 

The first two terms on the right hand side of (14) represent the organized periodic energy 
contribution, while the last term represents the turbulent energy contribution. 

Figure 2 shows the distribution of the mean axial velocity at the centerline. The anisotropic 
model gives better distribution downstream of the obstacle. Figure 3 compares the calculated 
distribution of the time-mean kinetic energy of the fluctuating motion (periodic + turbulent) 
along the centerline of the flow. The figure shows a better trend exhibited by the anisotropic 
k-£ model due to the improved resolution of the normal stresses. The standard k-£ model acts 
to damp the periodic fluctuations by producing too much eddy viscosity, which underestimates 
the time-averaged momentum transfer. Hence, the length of the separation region behind the 
obstacle is overpredicted. Also, the maximum of the kinetic energy at the centerline lies further 
downstream. The length of the recirculation zone and the location of the maximum fluctuating 
energy are improved by using the anisotropic model. The figure also shows some fluctuating 
energy in front of the obstacle, whereas measurements indicated that the flow remained virtually 
laminar there. This is because in the k-e model the large velocity gradients at the stagnation 
region produce large turbulent kinetic energy. Results are also obtained from calculations in 
which the production of k in front of the obstacle was suppressed. Figure 4 shows the mean 
axial velocity distribution indicating better comparison with the experiment downstream of the 
obstacle. Figure 5 shows the distribution of the mean kinetic energy along the centerline. It can 
be seen that suppressing the production of the kinetic energy in front of the obstacle causes an 
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0.16 0.10 



(a) Using the standard k-e model 
0.19 0.01 



(b) Using the anisotropic k-e model 
Fig. 6: < v* v* > /U 2 contours at T = 3 sec 

increase in the fluctuating energy. Also, the peak of the energy fluctuations is shifted slightly 
upstream closer to the experimental data. The figure also shows smaller residual fluctuating 
energy in front of the obstacle. Figure 6(a) and (b) show the contour plots of the normal 
turbulent stress term < v 1 v* > at an instant T = 3 sec. It can be seen that the anisotropic k-€ 
model produces higher < v* v 1 > values, which act to increase the total fluctuating energy. 

5. CONCLUSIONS 

The turbulent vortex shedding flow behind a square cylinder was analyzed using an efficient 
time-accurate numerical method based on the PISO methodology. Turbulence was modeled 
using an anisotropic k-e model which resolves the anisotropy of the Reynolds stresses rea- 
sonably well. Comparisons with the experimental data show the advantages of the model as 
compared with the standard isotropic k-e model. Accurate predictions, however, can only be 
made by accounting for the history and transport effects of the individual Reynolds stresses. 
The anisotropic k-e model seems to offer a compromise between the computationally intensive 
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Reynolds stress model and the standard isotropic k-e model. 
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Single point modeling of rotating turbulent flows 

By A. H. Hadid 1 , N. N. Mansour 2 AND O. Zeman 3 


A model for the effects of rotation on turbulence is proposed and tested. These 
effects which influence mainly the rate of turbulence decay are modeled in a modified 
turbulent energy dissipation rate equation that has explicit dependence on the mean 
rotation rate. An appropriate definition of the rotation rate derived from critical 
pQijyt; theory and b as ed on the invariants of the deformation tensor is proposed. 
The modeled dissipation rate equation is numerically well behaved and can be 
used in conjunction with any level of turbulence closure. The model is applied 
to the two-equation k-e turbulence model and is used to compute separated flows 
in a backward- facing step and an axisymmetric swirling coaxial jets into a sudden 
expansion. In general, the rotation modified dissipation rate model show some 
improvements over the standard k-e model. 

1. Motivation and objectives 

The ability to accurately model the effects of rotation on turbulence has a wide 
variety of important applications in rotating machinery and combustion devices. 
Many turbulent flows of engineering importance involve combinations of rotational 
and irrotational strains. However, turbulence models of the eddy viscosity type 
are oblivious to the presence of rotational strains since they depend only on the 
mean velocity gradients through their symmetric part (i.e. the mean rate of strain 
tensor). The rotation rate, for example, does not explicitly enter the equations for 
the turbulent kinetic energy and its dissipation rate, yet evidence from experiments 
(Wigeland and Nagib 1978, Jacquin et al. 1990) and from direct numerical simula- 
tion (Bardina et al. 1985, Speziale et al. 1987, Mansour et al. 1991) show that the 
decay rate of turbulence is reduced by the presence of rotation. 

The effects of rotation on turbulence are known to be subtle. They are manifested 
through changes in the spectr um of the turbulence caused by nonlinear interactions. 
For initially isotropic turbulence, rotation inhibits the cascade of energy from large 
to small scales. Zeman (1994) proposed a modified energy spectrum that takes into 
account the effects of rotati on at high Reynolds number by introducing a rotation 
wavenumber, kn = <y/f2 3 /e, below which rotation effects on spectral transfer are 
important. Much of the application work in simulating rotating flows have been 
conducted using varieties of eddy viscosity models (k-e or k-l) and second order 
closure models with modified dissipation rate transport equation to account for 
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rotational effects. However, most of these models fed to predict the asjmrptotrc 
of the turbulence decay rate in the hunts of large rotation rate. The 
objectives of this work are to model the effects of rotatton ustng single-point two 
“Son models and to offer an appropriate definition of the mean rotatton rate 
^consistent with the fact that spin is the mam cause of reductton m the 

dissipation rate. 


2. Accomplishments 

For incompressible viscous flow with constant properties, the 
eouations for the turbulent kinetic energy, t, and its dissipation rate, e, 
widely used for engineering applications take the form; 


Jb, t + Ujkj = D k + Pk-e 


£,. + U,,., =D, + P.-*. (2) 

where D„ and D, are the diffusion terms for k and e respectively and are modeled 




where „ is the laminar viscosity and i>, is the eddy viscosity = C„fc 2 /e. c k ««1 £ 
are the ratio of Prandtl to s chnndt jiui^rs and^e t Reynolds stress 

production term for k given as P k = wnere U < U J 

term and U, is the mean velocity in the i-direction. , 

16 Assuming that the production of th, nj , , P. . 

- "on rr - 

sipation rate is proportional to the turbulent energy dissipation rate term, i. . 
~ e /T. The modeled form of the dissipation rate equation becomes 


£ C 

t,t + Uje,j = D { + Cy-^Pk - C 2 -£ 


Due to the symmetry of the Reynolds stress tensor ujiij , the kinetic energy pro- 
duction term can be written as Pk = — UjUjSy, where S«j the* standard dis- 
the mean rate of strain tensor. Therefor it can be seen that 
sination rate eq. (3), has no explicit dependence on the mean rotation tensor 
q. = m - Uj i)/ 2. It follows that the commonly used modeled dissipa ion ra e 
^ition in only be affected indimctly by rotational strains through the chmtges 

that they induce in the Reynolds stress tensor. . j 

In order to sensitize the dissipation rate equation to rotational effects, consid^ 
the simple case of isotropic turbulence in a rotating frame. In t s case, an l 
decaying isotropic turbulence is described by; 


k, t = -« 
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e <* = ~ C *J ( 5 ) 

Equations (4) and (5) do not distiguish the difference between isotropic turbulence 
in a rotating frame and in an inertial frame. Models that have a non zero rotational 
correction have been proposed by Bardina ti al. (1985), for example, for rotating 
isotropic turbulence where eq. (5) takes the form 

e 2 

— C2 — Csfte (6) 


with C 2 = 1*83 and C 3 = 0.15. 

The above model is able only to accurately predict the reduction in the decay 
rate of the turbulent kinetic energy in rotating isotropic turbulence for weak to 
moderate rotation rates where the effects are small. However, for sufficiently high 
rotation rates and long enough time, the model drastically underpredicts the decay 
rate of the turbulent kinetic energy. 

Hanjallic and Launder (1980) proposed a model for which the c- transport equation 
in rotating isotropic turbulence takes the form 

e,t = -C 2 j-C 3 tfk (7) 

where Ci — 1.92 and C 3 = 0.27. 

This model predicts unphysical behavior of negative dissipation rate at high ro- 
tation rates, thus violating the realizability constraint. Other modifications to the 
dissipation rate transport equation have been proposed to account for rotational 
strains, e.g Raj (1975) and Pope (1978). Again they fail in one way or another to 
account accurately for the rotational effects. 

3. Proposed model 

In the present work a new model is proposed that accounts for rotational effects 
and correctly predicts the asymptotic behavior at zero to inifinte rotation rates. 
Consider the dissipation rate equation in rotating isotropic turbulence 


with 



a 2 \ e 2 
a 2 + 1 / k 


( 8 ) 


a = 0.35 Ro- 1 (9) 

where Ro is the Rossby number defined as Ro~ l = ftAr/e. For ft >> 1 , C 2 = 2.5, 
which gives a power law exponent n = 0.6 (in k ~ t ~ n ) maching the power law 
proposed by Squires et al. (1993) for the asymptotic state of rotating homogeneous 
turbulence at high Reynolds numbers. 

The experminental data of Jacquin et al. (1990) are used to test the proposed 
model. Their experiments consisted of measuring the velocity field and characteris- 
tic quantities characterizing the fluctuating field downstream of a rotating cylinder 


424 


A. E. Ha did, N. N. Mansour & O. Zeman 


£ 


containing a honycomb structure and a turbulence producing grid. The coupled 
differential equations for k and e describing the effects of rotation on an initially 
isotropic turbulence can be written as 

k ti = c (10) 

''- = -( C!+Cs ^t)t < u > 

These equations were solved numerically using a fourth-order Runge-Kutta inte- 
gration scheme. The model predictions (with C*2 = 1.7 and C$ = 5/6) are compared 
with the experimental data of Jacquin et al (1990) as shown in Fig. la. The model 
predicts well the evolution of turbulent kinetic energy and its decay rate for a wide 
range of rotation rates. We have also tested the model for the three Reynolds 
numbers measured by Jacquin et al (1990), and found similar agreement of the 
model predictions with the data. We should point out at this point that the value 
C*2 = 1-7, proposed here for zero rotation rate, is lower than the value convention- 
ally used in k-e modeling. We find that with the conventional value of C2 = 1.92 

(and C3 = 3/5) the model fails to predict the experimental data (see Fig. lb) 



FIGURE 1. Decay of turbulent kinetic energy. Symbols are the data of Jacquin 

et al (1990), lines are the model predictions, o h ft = 62.8 (rad/s), □ & 

ft = 31.4 (rad/s), * &; ft = 15.7. (a) Model predictions with C 2 = 1.7 

and C3 =5/6; (b) Model predictions with C 2 = 1.92 and C 3 = 3/5. 

4. Rotation Rate For General Flows 

In order to test the rotational correction proposed in eq. (8) to the dissipation 
rate equation for general flows where the rotation rate is a function of position and 
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in the presence of mean strains, the question arises as to what is the appropriate 
definition of the rotation rate, fi? 

In most previous studies, the rotation rate or the mean vorticity 0, was replaced 
by y/^jhij/2, where il.j = (U,j - U jt i)/2 is the rotation rate tensor of the mean 
flow. However, such definition does not distinguish between a vortex sheet and 
a vortex. A definition of a vortex or a region of vorticity (with spin) was given 
by Chong et al. (1990) -using the arguments of the critical point theory and the 
invariants of the deformation tensor- as a region in space where the vorticity is 
sufficiently strong to cause the rate of strain tensor to be dominated by the rotation 
tensor, i.e. the rate of deformation tensor has complex eigenvalues. This definition 
satisfies the principle of frame invariance since it depends only on the properties 
of the deformation tensor. We shall use it because the reduction in the dissipation 
rate is due mainly to the spin that the mean imposes on the turbulence. Consider 
the matrix D tJ of the elements of the deformation tensor, 

D tJ = U u (12) 

which can be split to 

Da = Si, -I- Uij ( 13 ) 

The complex eigenvalues of D tj axe found by solving the characteristic equation 
|£) tj _ A<5,j | = 0, where the A’s are the eigenvalues of Dij. For a 3 x 3 matrix, A can 
be found from the solution of 

A 3 + PA 2 4- Q\ + .R = 0 (14) 

where P, Q and R axe the matrix invariants and are given by 

P = -U it i ( 1 5 ) 
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then the three roots of A are; 


A + B, — 


A + B 
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+ i 


A-B 
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A- 


A + B 
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That is A can have: 

(:) all real roots which are distinct when 

[(<?/ 3) 3 +(H/2) 2 ] <0, 


or 

(it) ail real roots where at least two roots are equal when 


[(Q/3) 3 +(JJ/2) 2 ] =0, 


or 

(iit) one real root and a pair of complex conjugate roots when 

[(Q/3) 3 + (-R/2) 2 ] > 0. 


We shall follow Chong et a l. (1990) and define the rotation rate 

ft = 3(A) = ^(A - B ), when [(Q/3) 3 + (-R/2) 2 ] > 0, (19) 

ft = 0 otherwise. It is important to note that for two dimensional Cartesian flows, 
the rotation rate defined by Eq. (19) reduces to ft = y/]Q\, when Q, the determinant 
of the deformation tensor matrix, is negative. For pure shear the definition, eq. (19) 
yields ft = 0. Conventional models that are calibrated for shear flows, need not be 
recalibrated when corrections based on ft are added to the model. 


5. Numerical Procedure 

For a two-dimensional, incompressible and steady turbulent flow, the Reynolds 
averaged moment um , continuity, turbulent kinetic energy and dissipation rate equa- 
tions can be written in the generalized form; 


d_ 

dx 


( P U*) + i JVv*) - 1 (r,. § ) + (rr.. f ) + * (20) 


Where r = 1 for Cartesian two-dimensional flow, and y = r for two-dimensional ax- 
isymmetric flow. Table 1 gives a summary of the terms in eq. (20) for the dependent 
variables solved in the code. 
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Table 1. S umm ary of the governing equations, p is the density, T*. and r<t r 
are the exchange coefficients in the axial and radial directions respectively, S* is 
the source term for the variable $. In the table, p t is the effective viscosity given 
as p t = p + Pt, where p is the laminar viscosity and p t is the turbulent viscosity, 

Pt = C^pP/e. 

In the standard k-e turbulence model the constants C M , Ci, C 2 , <*k and a ( have 
the values 0.09, 1.44, 1.92, 1.0 and 1.0 respectively. 

In the rotation modified k-e turbulence model, only C 2 takes the form given by 

eq. (11) i.e, C 2 = 1-7 + (5/6 )a 2 /(a 2 + 1). 

The governing transport eq. (20) is solved using the primitive variables on a 
nonstaggered mesh and converted into a system of algebraic equations by integrat- 
ing over control volumes defined around each grid point. The SIMPLE pressure- 
correction scheme (Patankar 1980) is used to couple the pressure and velocities 
and the resulting algebraic equations are solved iteratively. The convective terms 
are differenced using a second-order upwind scheme while the diffusion terms are 
approximated by a central differencing scheme.The physical domain is discretized 
using a non-uniform mesh where grid points are clustered close to the walls. 

6. Model Application 

The performance of the present model for complicated recirculating flows is 
demonstrated through calculations and comparisons with the experimental data 
of Driver k. Seegmiller (1985) for backward-facing step flows and with the experi- 
ments of Roback k Johnson (1983) for a confined swirling coaxial jets into a sudden 

expansion. 

Figure 2, shows the streamlines for the backward-facing step using the rotation 
modified k-e turbulence model. The calculations were performed on a 100x4C) gn 
points. The computational domain had a length of 50# (H is the step height) an 
a width of 9 H. The experimental data were used to specify the inflow conditions 
for a channel flow calculation where the fully developed profiles at the channel exit 
were used as the inlet conditions for the backward-facing step calculations. Fully 
developed flow conditions were imposed at the outflow boundary. The stan ar 
wall function approach (Launder & Spalding 1974) was used to bridge the viscous 

sublayer near the wall. 
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FIGURE 2. Backward-facing step geometry and stream-function contours. The 
contour levels were set between (-0.1 and 0.1) with an increment level = 0.01. 
negative values, positive values. 

The computed reattachment lengths were 5.50# using the standard k-e turbu- 
lence model and 6.22# for the rotation modified k-e turbulence model. The modified 
k-e model prediction is closer to the experimental value of 6.10#. While these re- 
sults axe encouraging, they are mainly due to the fact that we have changed the 
value of C 2 for the non-rotating case. In general, a change in the value of C 2 will 
result in poor predictions of the mean profiles. The mean velocity profile at three 
locations downstream are shown on Fig. 3, while the turbulent stress profiles at 
X/H = 4 sire shown on Fig. 4. All the quantities were normalized by the step 
height (#) and the experimental reference free- stream velocity (?7 rc /). It can be 
seen that the overall performance of the rotation modified dissipation rate equation 
is better than the standard k-e model especially in the recirculation region (Figs. 3a, 
and 4). Some improvements are also obtained in the recovery region using the mod- 
ified k-e model. Figure 5 shows the contours of the effective rotation rate used as 
defined by Eq. (19). 

For the 2D/axisymmetric swirling flow computations, the expressions for the 
invariants Q and R (Eqs. (16) & (17) respectively) are expanded and Eq. (19) is 
used to obtain the values of fl. The model was used to predict the mean profiles for a 
confined double concentric jets with a swirling outer jet flow into a sudden expansion 
(Roback & Johnson, 1983, see Fig. 6). Measurements are available for the mean 
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FIGURE 3. Mean axial velocity profiles at different axial locations . o data (Driver 
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X/H = 4, (b) XI H = 8, (c) XI H = 12. 



0u03 004 *0020 -OJU13 Atuw 


(a) JR/cr?., (b) M <l v ~i 


FIGUEE 4. Tbrbulent stress profiles at X/H = 4. o data (Driver^Seegmiller 
19 g 5 ). modified k-e model; standard k-e model, (a) u l u 1 /U re} , (Dj 
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velocity profiles and velocity fluctuations downstream of the expansion. Simulations 
with a coarse nonuniform grid of 30x20 mesh points were made. However, there 
is some uncertainty about the inlet conditions to be used since the first velocity 
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FIGURE 5. Contours of the effective rotation rate, Cl. Contour levels were set 
between (0.1, 1.0) with an increment level = .01. are 



Vane S wirier 

Figure 6. Roback k Johnson’s swirling coaxial jets discharging into an expanded 
duct. 



profiles me as ured were 5 nun downstream of the expansion. 

To predict this flow, the measured profiles at 5mm were adjusted near the edges 
and were used as inlet conditions at the expansion plane. Preli min ary results ob- 
tained with the coarse mesh indicate similar trends as the experiment. Figure 7 
shows the streamline contours using the standard and the modified k-e turbulence 
models. The figure shows that a closed internal recirculation zone forms at the 
center with an additional zone at the comers downstream of the step. This causes a 
flow diversion outwards with high gradients between these regions. Figure 8 shows 
the axial and tangential velocity profiles at 25 mm downstream of the expansion 
using the standard and the modified k-e turbulence models. Results in this case 
indicate little or no improvements offered using the modified k-e model over the 
standard k-e model. Finer mesh may improve the results but the uncertainties in 
the inlet boundary conditions raise the question about the adequacy of using this 
experiment for validation purposes. 
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FIGURE 7. Swirling coaxial jets discharging into an expanded duct. Stream- 

function contour. levels were set between (-0.15,0.) with an increment level 

= 0 01j levels were set between (0.,0.7) with an increment level = 0.05. (a) 

Standard k-t model, (b) Modified k-e model. 



( a ) 



(b) 


FIGURE 8. Velocity profiles at X = 25 mm. o data (Roback t Jolmson 1983); 
modified k-e model; standard k-e model, (a) Axial Velocity, (b) T 

gential velocity. 

7. Conclusions 

A new simple model for the turbulent energy dissipation rate equation has been 
proposed to account for the rotational effects on turbulence. A frame invariant 
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definition of the rotation rate proposed by Chong et al. (1990) based on the critical 
point theory was used. The model can be used in conjunction with any level of 
| turbulence closure. It was applied to the two-equation k-e turbulence model and was 

tested for separted flows in a backward-facing step and for axisymmetric swirling jet 
into a sudden expansion. The model is numerically stable and showed improvements 
over the standard k-e turbulence model. It is important to point out that the 
present study was carried out to roughly evaluate the model, but that a systematic 
! recalibration of the constants in the k-e model is needed before going any further 

with the proposed model. 

The authors would like to acknowledge many discussions with Dr. K. Shariff 
j regarding proper definition of the rotation rate. 
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